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Abstract 

This document consists of lecture notes for a graduate course, which focuses on the 
relations between Information Theory and Statistical Physics. The course is aimed at 
EE graduate students in the area of Communications and Information Theory, as well 
as to graduate students in Physics who have basic background in Information Theory. 
Strong emphasis is given to the analogy and parallelism between Information Theory 
and Statistical Physics, as well as to the insights, the analysis tools and techniques that 
can be borrowed from Statistical Physics and 'imported' to certain problem areas in 
Information Theory. This is a research trend that has been very active in the last few 
decades, and the hope is that by exposing the student to the meeting points between 
these two disciplines, we will enhance his/her background and perspective to carry out 
research in the field. 

A short outline of the course is as follows: Introduction; Elementary Statistical 
Physics and its Relation to Information Theory; Analysis Tools in Statistical Physics; 
Systems of Interacting Particles and Phase Transitions; The Random Energy Model 
(REM) and Random Channel Coding; Additional Topics (optional). 
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1 Introduction 



This course is intended to EE graduate students in the field of Communications and Informa- 
tion Theory, and also to graduates of the Physics Department (in particular, graduates of the 
EE-Physics program) who have basic background in Information Theory, which is a prereq- 
uisite to this course. As its name suggests, this course focuses on relationships and interplay 
between Information Theory and Statistical Physics - a branch of physics that deals with 
many-particle systems using probabilitistic/statistical methods in the microscopic level. 

The relationships between Information Theory and Statistical Physics (-1- thermodynam- 
ics) are by no means new, and many researchers have been exploiting them for many years. 
Perhaps the first relation, or analogy, that crosses our minds is that in both fields, there 
is a fundamental notion of entropy. Actually, in Information Theory, the term entropy 
was coined after the thermodynamic entropy. The thermodynamic entropy was first intro- 
duced by Clausius (around 1850), whereas its probabihstic-statistical interpretation is due 
to Boltzmann (1872). It is virtually impossible to miss the functional resemblance between 
the two notions of entropy, and indeed it was recognized by Shannon and von Neumann. 
The well-known anecdote on this tells that von Neumann advised Shannon to adopt this 
term because it would provide him with "... a great edge in debates because nobody really 
knows what entropy is anyway. " 

But the relationships between the two fields go far beyond the fact that both share the 
notion of entropy. In fact, these relationships have many aspects, and we will not cover all 
of them in this course, but just to give the idea of their scope, we will mention just a few. 

• The Maximum Entropy (ME) Principle. This is perhaps the oldest concept that ties 
the two fields and it has attracted a great deal of attention, not only of information 
theortists, but also that of researchers in related fields like signal processing, image 
processing, and the like. It is about a philosopy, or a belief, which, in a nutshell, is the 
following: If in a certain problem, the observed data comes from an unknown probabil- 
ity distribution, but we do have some knowledge (that stems e.g., from measurements) 
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of certain moments of the underlying quantity /signal/random-variable, then assume 
that the unknown underlying probability distribution is the one with maximum entropy 
subject to (s.t.) moment constraints corresponding to this knowledge. For example, 
if we know the first and the second moment, then the ME distribution is Gaussian 
with matching first and second order moments. Indeed, the Gaussian model is perhaps 
the most widespread model for physical processes in Information Theory as well as 
in signal- and image processing. But why maximum entropy? The answer to this 
philosophical question is rooted in the second law of thermodynamics, which asserts 
that in an isolated system, the entropy cannot decrease, and hence, when the system 
reaches equilibrium, its entropy reaches its maximum. Of course, when it comes to 
problems in Information Theory and other related fields, this principle becomes quite 
heuristic, and so, one may question its relevance, but nevertheless, this approach has 
had an enormous impact on research trends throughout the last fifty years, after being 
proposed by Jaynes in the late fifties of the previous century, and further advocated 
by Shore and Johnson afterwards. In the book by Cover and Thomas, there is a very 
nice chapter on this, but we will not delve into this any further in this course. 

• Landauer's Erasure Principle. Another aspect of these relations has to do with a piece 
of theory whose underlying guiding principle is that information is a physical entity. In 
every information bit in the universe there is a certain amount of energy. Specifically, 
Landauer's erasure principle (from the early sixties of the previous century), which 
is based on a physical theory of information, asserts that every bit that one erases, 
increases the entropy of the universe by k In 2, where k is Boltzmann's constant. It is 
my personal opinion that these kind of theories should be taken with a grain of salt, 
but this is only my opinion. At any rate, this is not going to be included in the course 
either. 

• Large Deviations Theory as a Bridge Between Information Theory and Statistical Physics. 
Both Information Theory and Statistical Physics have an intimate relation to large de- 
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viations theory, a branch of probability theory which focuses on the assessment of the 
exponential rates of decay of probabilities of rare events, where the most fundamental 
mathematical tool is the Chernoff hound. This is a topic that will be covered in the 
course and quite soon. 

• Random Matrix Theory. How do the eigenvalues (or, more generally, the singular val- 
ues) of random matrices behave when these matrices have very large dimensions or if 
they result from products of many randomly selected matrices? This is a hot area in 
probability theory with many applications, both in Statistical Physics and in Infor- 
mation Theory, especially in modern theories of wireless communication (e.g., MIMO 
systems). This is again outside the scope of this course, but whoever is interested 
to 'taste' it, is invited to read the 2004 paper by Tulino and Verdii in Foundations 
and Trends in Communications and Information Theory, a relatively new journal for 
tutorial papers. 

• Spin Glasses and Coding Theory. It turns out that many problems in channel coding 
theory (and also to some extent, source coding theory) can be mapped almost ver- 
batim to parallel problems in the field of physics of spin glasses - amorphic magnetic 
materials with a high degree of disorder and very complicated physical behavior, which 
is cusomarily treated using statistical-mechanical approaches. It has been many years 
that researchers have made attempts to 'import' analysis techniques rooted in statis- 
tical physics of spin glasses and to apply them to analogous coding problems, with 
various degrees of success. This is one of main subjects of this course and we will 
study it extensively, at least from some aspects. 

We can go on and on with this list and add more items in the context of these very 
fascinating meeting points between Information Theory and Statistical Physics, but for now, 
we stop here. We just mention that the last item will form the main core of the course. We 
will see that, not only these relations between Information Theory and Statistical Physics 
are interesting academically on their own right, but moreover, they also prove useful and 
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beneficial in that they provide us with new insights and mathematical tools to deal with 
information-theoretic problems. These mathematical tools sometimes prove a lot more ef- 
ficient than traditional tools used in Information Theory, and they may give either simpler 
expressions for performance analsysis, or improved bounds, or both. 

At this point, let us have a brief review of the syllabus of this course, where as can be 
seen, the physics and the Information Theory subjects are interlaced with each other, rather 
than being given in two continuous, separate parts. This way, it is hoped that the relations 
between Information Theory and Statistical Physics will be seen more readily. The detailed 
structure of the remaining part of this course is as follows: 

1. Elementary Statistical Physics and its Relation to Information T/ieory; What is statis- 
tical physics? Basic postulates and the micro-canonical ensemble; the canonical en- 
semble: the Boltzmann-Gibbs law, the partition function, thermodynamical potentials 
and their relations to information measures; the equipartition theorem; generalized en- 
sembles (optional) ; Chernoff bounds and the Boltzmann-Gibbs law: rate functions in 
Information Theory and thermal equilibrium; physics of the Shannon limits. 

2. Analysis Tools in Statistical Physics: The Laplace method of integration; the saddle- 
point method; transform methods for counting and for representing non-analytic func- 
tions; examples; the replica method - overview. 

3. Systems of Interacting Particles and Phase Transitions: Models of many-particle sys- 
tems with interactions (general) and examples; a qualitative explanation for the ex- 
istence of phase transitions in physics and in information theory; ferromagnets and 
Ising models: the ID Ising model, the Curie- Weiss model; randomized spin-glass mod- 
els: annealed vs. quenched randomness, and their relevance to coded communication 
systems. 

4. The Random, Energy Model (REM) and Random Channel Coding'; Basic derivation and 
phase transitions - the glassy phase and the paramagnetic phase; random channel codes 
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and the REM: the posterior distribution as an instance of the Boltzmann distribution, 
analysis and phase diagrams, imphcations on code ensemble performance analysis. 

5. Additional Topics (optional): The REM in a magnetic field and joint source-channel 
coding; the generahzed REM (GREM) and hierarchical ensembles of codes; phase 
transitions in the rate-distortion function; Shannon capacity of infinite-range spin- 
glasses; relation between temperature, de Bruijn's identity, and Fisher information; 
the Gibbs inequality in Statistical Physics and its relation to the log-sum inequality 
of Information Theory. 

As already said, there are also plenty of additional subjects that fall under the umbrella 
of relations between Information Theory and Statistical Physics, which will not be covered 
in this course. One very hot topic is that of codes on graphs, iterative decoding, belief 
propagation, and density evolution. The main reason for not including these topics is that 
they are already covered in the course of Dr. Igal Sason: "Codes on graphs." 

I would like to emphasize that prior basic background in Information Theory will be 
assumed, therefore. Information Theory is a prerequisite for this course. As for the physics 
part, prior background in statistical mechanics could be helpful, but it is not compulsory. 
The course is intended to be self-contained as far as the physics background goes. The 
bibliographical list includes, in addition to a few well known books in Information Theory, 
also several very good books in elementary Statistical Physics, as well as two books on the 
relations between these two fields. 

As a final note, I feel compelled to clarify that the material of this course is by no means 
intended to be presented from a very comprehensive perspective and to consist of a full 
account of methods, problem areas and results. Like in many advanced graduate courses in 
our department, here too, the choice of topics, the approach, and the style strongly reflect 
the personal bias of the lecturer and his/her perspective on research interests in the field. 
This is also the reason that a considerable fraction of the topics and results that will be 
covered, are taken from articles in which I have been involved. 
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2 Elementary Stat. Physics and Its Relation to IT 

2.1 What is Statistical Physics? 

Statistical physics is a branch in Physics which deals with systems with a huge number 
of particles (or any other elementary units), e.g., of the order of magnitude of Avogadro's 
number, that is, about 10^^ particles. Evidently, when it comes to systems with such an 
enormously large number of particles, there is no hope to keep track of the physical state 
(e.g., position and momentum) of each and every individual particle by means of the classical 
methods in physics, that is, by solving a gigantic system of differential equations pertaining to 
Newton's laws for all particles. Moreover, even if these differential equations could have been 
solved (at least approximately), the information that they would give us would be virtually 
useless. What we normally really want to know about our physical system boils down 
to a bunch of macroscopic parameters, such as energy, heat, pressure, temperature, volume, 
magnetization, and the like. In other words, while we continue to believe in the good old laws 
of physics that we have known for some time, even the classical ones, we no longer use them in 
the ordinary way that we are familar with from elementary physics courses. Rather, we think 
of the state of the system, at any given moment, as a realization of a certain probabilistic 
ensemble. This is to say that we approach the problem from a probabilistic (or a statistical) 
point of view. The beauty of statistical physics is that it derives the macroscopic theory of 
thermodynamics (i.e., the relationships between thermodynamical potentials, temperature, 
pressure, etc.) as ensemble averages that stem from this probabilistic microscopic theory - 
the theory of statistical physics, in the limit of an infinite number of particles, that is, the 
thermodynamic limit. As we shall see throughout this course, this thermodynamic limit is 
parallel to the asymptotic regimes that we are used to in Information Theory, most notably, 
the one pertaining to a certain 'block length' that goes to infinity. 
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2.2 Basic Postulates and the Microcanonical Ensemble 



For the sake of concreteness, let us consider the example where our many-particle system is a 
gas, namely, a system with a very large number n of mobile particles, which are free to move 
in a given volume. The microscopic state (or microstate, for short) of the system, at each 
time instant t, consists, in this example, of the position fi{t) and the momentum Pi{t) of each 
and every particle, 1 < i < n. Since each one of these is a vector of three components, the 
microstate is then given by a (6n)-dimensional vector x(t) = {{fi(t),'Pi(t)), i = 1,2, . . . , n}, 
whose trajectory along the time axis, in the phase space, JR^^, is called the phase trajectory. 

Let us assume that the system is closed, i.e., isolated from its environment, in the sense 
that no energy flows inside or out. Imagine that the phase space IR^*^ is partitioned into very 
small hypercubes (or cells) Ap x Ar. One of the basic postulates of statistical mechanics is 
the following: In the very long range, the relative amount of time at which x{t) spends at 
each such cell converges to a certain number between and 1, which can be given the meaning 
of the probability of this cell. Thus, there is an underlying assumption of equivalence between 
temporal averages and ensemble averages, namely, this is the assumption of ergodicity. 

What are then the probabilities of these cells? We would like to derive these probabilities 
from first principles, based on as few as possible basic postulates. Our first such postulate 
is that for an isolated system (i.e., whose energy is fixed) all microscopic states {x{t)} are 
equiprobable. The rationale behind this postulate is twofold: 

• In the absence of additional information, there is no apparent reason that certain 
regions in phase space would have preference relative to any others. 

• This postulate is in harmony with a basic result in kinetic theory of gases - the Liouville 

theorem, which we will not touch upon in this course, but in a nutshell, it asserts that 

the phase trajectories must lie along hypersurfaces of constant probability densityj^ 

^This is a result of the energy conservation law along with the fact that probability mass behaves like 
an incompressible fluid in the sense that whatever mass that flows into a certain region from some direction 
must be equal to the outgoing flow from some other direction. This is reflected in the so called continuity 
equation. 
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Before we proceed, let us slightly broaden the scope of our discussion. In a more general 
context, associated with our n-particle physical system, is a certain instantaneous microstate, 
generically denoted by a; = (xi, X2, . . . , where each Xi, 1 < i < n, may itself be a vector 
of several physical quantities associated particle number i, e.g., its position, momentum, 
angular momentum, magnetic moment, spin, and so on, depending on the type and the 
nature of the physical system. For each possible value of x, there is a certain Hamiltonian 
(i.e., energy function) that assigns to a; a certain energy £^(a;),^ Now, let us denote by Vl{E) 
the density-of-states function, i.e., the volume of the shell {x : S{x) = E}, or, slightly 
more precisely, Q{E)dE = Vol{a; : E < S{x) < E + dE}, which will be denoted also as 
Vol{a; : S{x) ^ E}, where the dependence on dE will normally be ignored since Q{E) is 
typically exponential in n and dE will have virtually no effect on its exponential order as 
long as it is small. Then, our above postulate concerning the ensemble of an isolated system, 
which is called the microcanonincal ensemble, is that the probability density P{x) is given 



In the discrete case, things are, of course, a lot easier: Then, Q{E) would be the number 
of microstates with £{x) = E (exactly) and P{x) would be the uniform probability mass 
function across this set of states. In this case, Q{E) is analogous to the size of a type class 
in Information Theory, and P{x) is the uniform distribution across this type class. 

Back to the continuous case, note that Q{E) is, in general, not dimensionless: In the 
above example of a gas, it has the physical units of [length x momentum]^", but we must get 
rid of these physical units because very soon we are going to apply non-linear functions on 
^1{E), like the logarithmic function. Thus, we must normalize this volume by an elementary 
reference volume. In the gas example, this reference volume is taken to be h^"', where h 
is Planck's constant ^ 6.62 x 10^^^ Joules-sec. Informally, the intuition comes from the 
fact that h is our best available "resolution" in the plane spanned by each component of 

^For example, in the case of an ideal gas, £{x) = X)"=i ^^2m ' independently of the positions {ri}, namely, 
it accounts for the contribution of the kinetic energies only. In more complicated situations, there might be 
additional contributions of potential energy, which depend on the positions. 



by 





elsewhere 



(1) 
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fi and the corresponding component of Pi, owing to the uncertainty principle in quantum 
mechanics, which tells us that the product of the standard deviations Apa ■ Ar^ of each 
component a {a = x,y,z) is lower bounded by h/2, where h = /i/(27r). More formally, this 
reference volume is obtained in a natural manner from quantum statistical mechanics: by 
changing the integration variable ptokhj using p = hk, where k is the wave vector. This is 
a well-known relationship pertaining to particle-wave duality. Now, having redefined Q{E) 
in units of this reference volume, which makes it then a dimensionless quantity, the entropy 
is defined as 

S{E) = klnQ{E), (2) 

where k is Boltzmann's constant ^ 1.38 x 10^^^ Joule/degree. We will soon see what is the 
relationship between S{E) and the information-theoretic entropy. 

To get some feeling of this, it should be noted that normally, Q{E) behaves as an exponen- 
tial function of n (at least asymptotically), and so, S{E) is roughly linear in n. For example, 
if £{x) = X]r=i ^^1^' then ^l(E) is the volume of a shell or surface of a (3n)-dimensional 
sphere with radius y/^mE, which is proportional to {2mE)^"'^'^V"', but we should divide this 
by n\ to account for the fact that the particles are indistinguishable and we don't count 
permutations as distinct physical states in this casej^l More precisely, one obtains: 



S{E) = k\n 



f A-nmEV'^''^ V 



I 1 • —7-:r- + -n/c n/c In I 1 

V 3n y nW"" 2 V 3n y 



fATTmEV^ V 



+ \nk. (3) 



Assuming E oc n and V oc n, we get S{E) (x n. A physical quantity like this, that has 
a linear scaling with the size of the system n, is called an extensive quantity. So, energy, 
volume and entropy are extensive quantities. Other quantities, which are not extensive, i.e., 
independent of the system size, like temperature and pressure, are called intensive. 

It is interesting to point out that from the function S{E), or actually, the function 
S{E,V,n), one can obtain the entire information about the relevant macroscopic physical 



•^Since the particles are mobile and since they have no colors and no identity certficiates, there is no 
distinction between a state where particle no. 15 has position r and momentum p while particle no. 437 has 
position r' and momentum jf and a state where these two particles are swapped. 
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quantities of the system, e.g., temperature, pressure, and so on. The temperature T of the 
system is defined according to: 

1 fdS{E)\ 



T 



V dE ) 



V 



(4) 

where {■)v means that the derivative is taken in constant volume^ Intuitively, in most 
situations, we expect that S{E) would be an increasing function of E (although this is not 
strictly always the case), which means T > 0. But T is also expected to be increasing with 
E (or equivalently, E is increasing with T, as otherwise, the heat capacity dE/dT < 0). 
Thus, 1/T should decrease with E, which means that the increase of S" in slows down 
as E grows. In other words, we expect S{E) to be a concave function of E. In the above 
example, indeed, S{E) is logarithmic in E and we get 1/T = dS/dE = 3nk/{2E), which 
means E = 3nkT/2. Pressure is obtained by P = T ■ dS/dV, which in our example, gives 
rise to the state equation of the ideal gas, P = nkT/V. 

How can we also see mathematically that under "conceivable conditions", S{E) is a 
concave function? We know that the Shannon entropy is also a concave functional of the 
probability distribution. Is this related? 

As both E and S are extensive quantities, let us define E = ne and 

s(e) = hm (5) 

i.e., the per-particle entropy as a function of the per-particle energy. Consider the case 
where the Hamiltonian is additive, i.e., 

n 

Six) = J2n^^) (6) 

i=l 

just like in the above example where £{x) = ^1^=1 ''^l^- Then, obviously, 

Q{niei + 72262) > ri(niei) ■ (^(^262), (7) 



This definition of temperature is related to tfie classical thcrmodynamical definition of entropy as 
dS = dQ/T, where Q is heat, as in the absence of external work, when the volume V is fixed, all the energy 
comes from heat and so, dE = dQ. 
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and so, we get: 



k \nQ{niei + ^262) 



> 



klnQ{niei) A; In (77-2 62) 



rii + n2 ni+ n2 

ni In r2(niei) 



k In VL{n2e2) 



(8) 



ni + n2 rii 



rii + n2 



and so, by taking rii and n2 to 00, with rii/ {rii + 722) — ?■ A G (0, 1), we get: 



s(Aei + (1 - A)e2) > As(ei) + (1 - A)s(e2), 



(9) 



which estabhshes the concavity of s(-) at least in the case of an additive Hamiltonian, which 
means that the entropy of mixing two systems of particles is greater than the total entropy 
before they are mixed (the second law). A similar proof can be generalized to the case 
where S{x) includes also a limited degree of interactions (short range interactions), e.g., 
S{x) = J2'i=i^i^ij^i+i)j but this requires somewhat more caution. In general, however, 
concavity may no longer hold when there are long range interactions, e.g., where some 
terms of S{x) depend on a linear subset of particles. Simple examples can be found in: 
H. Touchette, "Methods for calculating nonconcave entropies," arXiv:1003.0382vl [cond- 
mat.stat-mech] 1 Mar 2010. 

Example - Schottky defects. In a certain crystal, the atoms are located in a lattice, and at 
any positive temperature there may be defects, where some of the atoms are dislocated (see 
Fig. [1]). Assuming that defects are sparse enough, such that around each dislocated atom all 
neighors are in place, the activation energy, eo, required for dislocation is fixed. Denoting 
the total number of atoms by and the number of defected ones by n, the total energy is 
then E = neo, and so. 




(10) 



or, equivalently. 



S{E) 




k[N In A^ — n Inn — (A^ — n) ln{N — n)] by the Stirling approximation 
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Thus, 

f ~ dE ~ 'dm'dE ~ 7^' n ' ^ ^ 

which gives the number of defects as 

exp(eo/A;r) + l- ^^^^ 
At T = 0, there are no defects, but their number increases gradually with T, approximately 
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Figure 1: Schottky defects in a crystal lattice. 

according to exp(— eo/ZcT). Note that from a shghly more information-theoretic point of 
view, 

S{E) = Hn ( ^ ) kNh, (^) = kNh, (^) = kNh, (l) , (13) 



where 



h2{x) = —xlnx — {1 — x) ln(l — x). 



Thus, the thermodynamical entropy is intimately related to the Shannon entropy. We will 
see shortly that this is no coincidence. Note also that S{E) is indeed concave in this example. 
□ 

What happens if we have two independent systems with total energy E, which lie in 
equilibrium with each other. What is the temperature T? How does the energy split between 
them? The number of combined microstates where system no. 1 has energy Ei and system 
no. 2 has energy E2 — E — Ei is Qi{Ei) ■ Q2{E — Ei). If the combined system is isolated, 
then the probability of such a combined microstate is proportional to Jli(-E'i) ■Q,2{E — Ei). 
Keeping in mind that normally, Vti and Vt2 are exponential in n, then for large n, this 
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product is dominated by the value of Ei for which it is maximum, or equivalently, the sum 
of logarithms, Si{Ei) + S2{E — Ei), is maximum, i.e., it is a mctximum entropy situation, 
which is the second law of thermodynamics. This maximum is normally achieved at 
the value of Ei for which the derivative vanishes, i.e., 

S[{E,) - S',{E - E,) = (14) 

or 

S[{E,)-S',{E2) = Q (15) 

which means 

^ ^ S[{E,) = S',{E,) ^ 1. (16) 

J-l J-2 

Thus, in equilibrium, which is the maximum entropy situation, the energy splits in a way 
that temperatures are the same. 

2.3 The Canonical Ensemble 

So far we have assumed that our system is isolated, and therefore has a strictly fixed energy 
E. Let us now relax this assumption and assume that our system is free to exchange energy 
with its large environment (heat bath) and that the total energy of the heat bath Eq is by 
far larger than the typical energy of the system. The combined system, composed of our 
original system plus the heat bath, is now an isolated system at temperature T. So what 
happens now? 

Similarly as before, since the combined system is isolated, it is governed by the micro- 
canonical ensemble. The only difference is that now we assume that one of the systems (the 
heat bath) is very large compared to the other (our test system). This means that if our 
small system is in microstate x (for whatever definition of the microstate vector) with energy 
£{x)^ then the heat bath must have energy Eq — E{x) to complement the total energy to 
Eq. The number of ways that the heat bath may have energy Eq — 8{x) is Qhb{Eo — ^i^)), 
where 0//b(-) is the density-of-states function pertaining to the heat bath. In other words. 
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the number of microstates of the combined system for which the small subsystem is in mi- 
crostate x is Qhb{Eq — S{x)). Since the combined system is governed by the microcanonical 
ensemble, the probability of this is proportional to VIhb{,Eq — £{x)). More precisely: 

nHB{Eo-£{x)) 



^^^^ Y.^,nHB{E^-s{x')y 

Let us focus on the numerator for now, and normalize the result at the end. Then, 

P{x) oc VLhb{Eq-S{x)) 

= e-K^{SHB{E^-S{x))/k} 

Shb{Eo) ldSHB{E) 



(17) 



exp 



exp 



k 

Shb{Eo) 
k 



k 
1 



dE 
Six) 



8{x) 



E=Eq 



oc exp{-£:(a;)/(A;T)}. 



It is customary to work with the so called inverse temperature: 



and so. 



^-kf 



P(x) oc e-'^^(^). 



(18) 



(19) 



(20) 



Thus, all that remains to do is to normalize, and we then obtain the Boltzmann-Gibbs (B-G) 
distribution, or the canonical ensemble, which describes the underlying probability law in 
equilibrium: 



p/ X _ e^p{-l3£{X)} 



where Z[I3) is the normalization factor: 

Z(^)=5]exp{-^£:(a;)} 

in the discrete case, or 



X 



Z{/3)^ J dxexp{-/3£{x)} 



(21) 



(22) 



18 



in the continuous case. 

This is one of the most fundamental results in statistical mechanics, which was obtained 
solely from the energy conservation law and the postulate that in an isolated system the 
distribution is uniform. The function Z{f3) is called the partition function, and as we shall 
see, its meaning is by far deeper than just being a normalization constant. Interestingly, a 
great deal of the macroscopic physical quantities, like the internal energy, the free energy, the 
entropy, the heat capacity, the pressure, etc., can be obtained from the partition function. 

The B-G distribution tells us then that the system "prefers" to visit its low energy states 
more than the high energy states. And what counts is only energy differences, not absolute 
energies: If we add to all states a fixed amount of energy Eq, this will result in an extra 
factor of e~^^° both in the numerator and in the denominator of the B-G distribution, which 
will, of course, cancel out. Another obvious observation is that whenever the Hamiltonian 
is additive, that is, S{x) — Yli=i^{^i)j the various particles are statistically independent: 
Additive Hamiltonians correspond to non-interacting particles. In other words, the {xi}'s 
behave as if they were drawn from a memoryless source. And so, by the law of large numbers 
^Yl^i=i^ixi) will tend (almost surely) to e = E{£{Xi)}. Nonetheless, this is different from 
the microcanonical ensemble where ^Y17=i^i-^i) held strictly at the value of e. The 
parallehsm to Information Theory is as follows: The microcanonical ensemble is parallel 
to the uniform distribution over a type class and the canonical ensemble is parallel to a 
memoryless source. 

The two ensembles are asymptotically equivalent as far as expectations go. They continue 
to be such even in cases of interactions, as long as these are short range. It is instructive 
to point out that the B-G distribution could have been obtained also in a different manner, 
owing to the maximum-entropy principle that we mentioned in the Introduction. Specifically, 
consider the following optimization problem: 

max H{X) 

s.t. '^P{x)S{x) = E [or in physicists' notation: {S{X)) = E] (23) 

X 
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By formalizing the equivalent Lagrange problem, where ^ now plays the role of a Lagrange 
multiplier: 



max {h{X) + I3 



X 



(24) 



or equivalently, 



min |5]P(a.)£(a.)-^| (25) 

one readily verifies that the solution to this problem is the B-G distribution where the 
choice of /3 controls the average energy E. In many physical systems, the Hamiltonian is 
a quadratic (or "harmonic") function, e.g., |mv^, ^kx"^, ^CV'^, ^LP, ^luj'^, etc., in which 
case the resulting B-G distribution turns out to be Gaussian. This is at least part of the 
explanation why the Gaussian distribution is so frequently encountered in Nature. Note 
also that indeed, we have already seen in the Information Theory course that the Gaussian 
density maximizes the (differential) entropy s.t. a second order moment constraint, which is 
equivalent to our average energy constraint. 

2.4 Properties of the Partition Function and the Free Energy 

Let us now examine more closely the partition function and make a few observations about 
its basic properties. For simphcity, we shall assume that x is discrete. First, let's look at 
the hmits: Obviously, Z(0) is equal to the size of the entire set of microstates, which is also 

^^0(£'), This is the high temperature limit, where all microstates are equiprobable. At 
the other extreme, we have: 

hm = - min S{x) ^ -Eos (26) 

/3-)-oo p X 

which describes the situation where the system is frozen to the absolute zero. Only states 
with minimum energy - the ground-state energy, prevail. 

Another important property of Z{P), or more precisely, of InZ(^), is that it is a log- 
moment generating function: By taking derivatives of InZ(^), we can obtain moments (or 
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cumulants) of S{X). For the first moment, we have 



E{S{X)} = {S{X)) - ^^^,,six) (27) 

Similarly, it is easy to show (exercise) that 

Var{^(X)} = {S\X)) - {S{X)f = (28) 

_|2 

This in turn implies that — ^ 0' which means that In Z{P) must always be a convex 
function. Higher order derivatives provide higher order moments. 

Next, we look at Z shghtly differently than before. Instead of summing e"'^^^^) across 
all states, we go by energy levels (similarly as in the method of types). This amounts to: 



X 



~ ^ gn.(6)/fc . g-/3n6 recall that 5(ne) fti ns(e) 

e 

= 5^exp{-n/3[e-Ts(e)]} 

= maxexp{— n^[e — Ts(e)]} 

= exp{—n/3 min[e — Ts(e)]} 

= exp{-n/3[e* -rs(e*)]} 

^ e-^^ (29) 



A 



The quantity / = e — ^^(e) is the (per-particle) free energy. Similarly, the entire free energy, 
F, is defined as 

F^E-TS^-^^. (30) 

The physical meaning of the free energy is this: A change, or a difference, AF = F2 — Fi, 
in the free energy means the minimum amount of work it takes to transfer the system 
from equilibrium state 1 to another equilibrium state 2 in an isothermal (fixed temperature) 
process. And this minimum is achieved when the process is quasistatic, i.e., so slow that 
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the system is always almost in equilibrium. Equivalently, —AF is the maximum amount of 
work that that can be exploited from the system, namely, the part of the energy that is free 
for doing work (i.e., not dissipated as heat) in fixed temperature. Again, this maximum is 
attained by a quasistatic process. 

We see that the value e* of e that minimizes /, dominates the partition function and 
hence captures most of the probability. As n grows without bound, the energy probability 
distribution becomes sharper and sharper around ne*. Thus, we see that equilibrium in the 
canonical ensemble amounts to minimum free energy. This extends the second law of 
thermodynamics from the microcanonical ensemble of isolated systems, whose equilibrium 
obeys the maximum entropy principle. The maximum entropy principle is replaced, more 
generally, by the minimum free energy principle. Note that the Lagrange minimization 
problem that we formalized before, i.e.. 



is nothing but minimization of the free energy, provided that we identify H with the physical 
entropy S (to be done very soon) and the Lagrange multiplier with kT. Thus, the B-G 
distribution minimizes the free energy for a given temperature. 

Although we have not yet seen this explicitly, but there were already hints and terminol- 
ogy suggests that the thermodynamical entropy S{E) is intimately related to the Shannon 
entropy H{X). Wc will also see it shortly in a more formal manner. But what is the 
information-theoretic analogue of the free energy? 

Here is a preliminary guess based on a very rough consideration: The last chain of 
equalities reminds us what happens when we sum over probabilities typc-by-type in IT 
problems: The exponentials exp{— /3£(a;)} are analoguous (up to a normalization factor) to 
probabilities, which in the memoryless case, are given by P{x) — exp{— + D{P\\P)]}. 
Each such probability is weighted by the size of the type class, which as is known from the 
method of types, is exponentially e"^, whose physical analogue is Q,{E) — e"'*^^)/*^. The 
product gives exp{— nL'(P||P)} in IT and exp{— n/?/} in statistical physics. This suggests 




(31) 
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that perhaps the free energy has some analogy with the divergence. Is this true? We will 
see shortly a somewhat more rigorous argument. 

More formally, let us define 

m = lim (32) 

n— ^cxD n 

and, in order to avoid dragging the constant k, let us define S(e) = lim„^oo ^ In ^^(^e) = 
s{e)/k. Then, the above chain of equalities, written slighlty differently, gives 



n— s>oo n 

,[S(e)-/3e] 



= lim - In J y e"[' 

= max[S(e) — /3e]. 

Thus, is (a certain variant of) the Legendre transforrrj^ oi S(e). As S(e) is (normally) 
a concave function, then it can readily be shown (execrise) that the inverse transform is: 

S(e) = min[/3e + 0(/3)]. (33) 

The achiever, e*(/3), of 0(/3) in the forward transform is obtained by equating the derivative 
to zero, i.e., it is the solution to the equation 

/3 = S'(e), (34) 

or in other words, the inverse function of S'(-). By the same token, the achiever, /3*(e), of 
S(e) in the backward transform is obtained by equating the other derivative to zero, i.e., it 
is the solution to the equation 

e = -</>'(/?) (35) 
or in other words, the inverse function of — 0'(-). 

Exercise: Show that the functions S'(-) and — </>'(■) are inverses of one another. □ 

This establishes a relationship between the typical per-particle energy e and the inverse 



^More precisely, the ID Legendre transform of a real function f{x) is defined as g{y) — snp^lxy — f{x)]. 
If / is convex, it can readily be shown that: (i) The inverse transform has the very same form, i.e., f{x) = 
swpy[xy — g{y)], and (ii) The derivatives f'{x) and g'{y) are inverses of each other. 
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temperature ^ that gives rise to e (cf. the Lagrange interpretation above, where we said that 
^ controls the average energy) . Now, obersve that whenever ^ and e are related as explained 
above, we have: 

E(6) = /3e + m = m - P ■ (36) 

On the other hand, if we look at the Shannon entropy pertaining to the B-G distribution, 

we get: 

H{X) = lim-£;|ln ^^^^ 



lim 

n— >oo 



lnZ(/3) ^ PE{£iX)} 



n n 
= </.(/?) -/?• </.'(/?). 

which is exactly the same expression as before, and so, E(e) and H are identical whenever 
P and e are related accordingly. The former, as we recall, we defined as the normalized 
logarithm of the number of microstates with per-particle energy e. Thus, we have learned 
that the number of such microstates is exponentially e**^, a result that looks familar to 
what we learned from the method of types in IT, using combinatorial arguments for finite- 
alphabet sequences. Here we got the same result from substantially different considerations, 
which are applicable in situations far more general than those of finite alphabets (continuous 
alphabets included). Another look at this relation is the following: 

^ exp{-/3ne-n(l){l3)} = n{ne) ■exp{-n[/3e + (f){/3)]} (37) 

X: £{X)Kine 

which means that VL{ne) < exp{n[/3e + 0(/9)]} for all /3, and so, 

fi(ne) < exp{nmm[^e + 0(^)]} = e""^^'^ = e"^. (38) 

A compatible lower bound is obtained by observing that the minimizing ^ gives rise to 
{£{Xi)) — e, which makes the event {x : S{x) fa ne} a high-probabihty event, by the 
weak law of large numbers. A good reference for further study and from a more general 
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perspective is: 

M. J. W. Hall, "Universal geometric approach to uncertainty, entropy, and information," 
Phys. Rev. A, vol. 59, no. 4, pp. 2602-2615, April 1999. 

Having established the identity between the Shannon-theoretic entropy and the thermo- 
dynamical entropy, we now move on, as promised, to the free energy and seek its information- 
theoretic counterpart. More precisely, we will look at the difference between the free energies 
of two different probability distributions, one of which is the B-G distibution. Consider first, 
the following chain of equalities concerning the B-G distribution: 



P{x) = 



ex.p{-\nZ{/3)- /3S{x)} 



= eMmP)-S{x)]}. (39) 

Consider next another probability distribution Q, different in general from P and hence 
corresponding to non-equilibrium. Let us now look at the divergence: 

Q{x) 



D{Q\\P) = J] g(a;) In 



X 



P{x) 



X 

= -HQ-f5Y,Q{^)[Fp-£{x)] 

X 



or equivalently. 



FQ^Fp + kT-D{Q\\P) 



Thus, the free energy difference is indeed related to the the divergence. For a given tem- 
perature, the free energy away from equilibrium is always larger than the free energy at 
equilibrium. Since the system "wants" to minimize the free energy, it eventually converges 
to the B-G distribution. More details on this can be found in: 
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1. H. Qian, "Relative entropy: free energy Phys. Rev. E, vol. 63, 042103, 2001. 

2. G. B. Bagci, arXiv:cond-mat/070300vl, 1 Mar. 2007. 

Another interesting relation between the divergence and physical quantities is that the di- 
vergence is proportional to the dissipated work (^average work — free energy difference) 
between two equilibrium states at the same temperature but corresponding to two different 

values of some external control parameter. Details can be found in: R. Kawai, J. M. R. Par- 
rondo, and C. Van den Broeck, "Dissipation: the phase-space perspective," Phys. Rev. Lett., 
vol. 98, 080602, 2007. 

Let us now summarize the main properties of the partition function that we have seen 
thus far: 

1. Z{/3) is a continuous function. Z{0) = and lim^^oo = —Egs- 

2. Generating moments: (S) = -dlnZ/d^, Var{£:(X)} = d^lnZ/dp^ convexity of 
InZ, and hence also of 0(/3). 

3. and E are a Legendre-transform pair. E is concave. 

4. E(e) coincides with the Shannon entropy of the B-G distribution. 

5. Fq = Fp + kT-D{Q\\P). 

Exercise: Consider Z{P) for an imaginary temperature (5 = ju, where j = a/— T, and define 
z{E) as the inverse Fourier transform of Z{juj). Show that z{E) — fl{E) is the density of 
states, i.e., for Ei < E2, the number of states with energy between Ei and E2 is given by 

jll z{E)dE. □ 

Thus, Z{-) can be related to energy enumeration in two different ways: one is by the Lcgcndrc 
transform of InZ for real /3, and the other is by the inverse Fourier transform of Z for 
imaginary (3. This double connection between Z and fl is no coincidence, as we shall see 
later on. 
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Example - A two level system. Similarly to the earlier example of Schottky defets, which 
was previously given in the context of the microcanonical ensemble, consider now a system 

of n independent particles, each having two possible states: state of zero energy and state 
1, whose energy is eo, i.e., S{x) = cqx, x e {0, 1}. The XiS are independent, each having a 
marginal: 



P{x) 



X e {0,1}. 



In this case. 



and 



1 + e-^«o 
0(^) = ln(l + e-^^° 



S(e) = min[/3e + ln(l + e-'^'o) 



To find P*{e), we take the derivative and equate to zero: 

-/3eo 



e 



I) 



which gives 



1 + e-'^^o 
ln(e/eo - 1) 



eo 



On substituting this back into the above expression of S(e), we get: 



E(e) = - In 
Co 



- 1 + In 



1 + exp < — In 



- 1 



eo 



which after a short algebraic manipulation, becomes 

E(e) = h2 

just like in the Schottky example. In the other direction: 



(f){f3) — max 



I3e 



whose achiever e*{f3) solves the zero-derivative equation: 



1 

eo 



In 



" 1 - e/ep " 
e/eo 



(40) 



(41) 



(42) 



(43) 



(44) 



(45) 



(46) 



(47) 



(48) 
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or equivalently, 

which is exactly the inverse function of f3*{e) above, and which when plugged back into the 
expression of (f){f3), indeed gives 

=ln(l + e-^'°). □ (50) 



Comment: A very similar model (and hence with similar results) pertains to non-interacting 
spins (magnetic moments), where the only difference is that x G { — 1,+1} rather than 
X G {0, 1}. Here, the meaning of the parameter eg becomes that of a magnetic field, which 
is more customarily denoted by B (or H) , and which is either parallel or antiparallel to that 
of the spin, and so the potential energy (in the appropriate physical units), B ■ x, is either 
Bx or —Bx. Thus, 

l3Bx 

^W = 2^^*(M^ Z(/3)=2cosh(/3B). (51) 
The net magnetization per-spin is defined as 



A 

m = 



(^|:-^-)=(^.> = gi^ = ta^(W (62) 



This is the paramagnetic characteristic of the magnetization as a function of the magnetic 
field: As B ^ ±00, the magnetization m — > ±1 accordingly When the magnetic field is 
removed {B — 0), the magnetization vanishes too. We will get back to this model and its 
extensions in the sequel. □ 

Exercise: Consider a system of n non-interacting particles, each having a quadratic Hamil- 
tonian, £{x) — ^ax^, x Show that here, 

E,)4.n(i^) ,53) 



and 



W)4ln(^). (54) 



Show that /3*(e) = l/(2e) and hence e*(/3) = 1/(2/3). 
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2.5 The Energy Equipartition Theorem 

Prom the last exercise, we have learned that for a quadratic Hamiltonian, £{x) — ^ax"^, we 
have €*{/3), namely, the average per-particle energy, is given l/{2/3) — kT/2, independently 
of a. If we have n such quadratic terms, then of course, we end up with nkT/2. In the 
case of the ideal gas, we have 3 such terms (one for each dimension) per particle, thus a 
total of 3n terms, and so, E — 3nkT/2, which is exactly what we obtained also in the 
microcanonical ensemble, which is equivalent (recall that this was obtained then by equating 
1/T to the derivative of S{E) = k ln[const x E^""/^]). In fact, we observe that in the canonical 
ensemble, whenever we have an Hamiltonian of the form ^xj+ some arbitrary terms that do 
not depend on then Xi is Gaussian (with variance kT/a) and independent of the other 
guys, i.e., p{xi) oc e""^^/^^*^^^. Hence it contributes an amount of 

to the total average energy, independently of a. It is more precise to refer to this Xi as a 
degree of freedom rather than a particle. This is because in the 3D world, the kinetic energy, 
for example, is given by p^/(2m) +p^/(2m) +p^/(2m), that is, each particle contributes 
three additive quadratic terms rather than one (just like three independent one-dimensional 
particles) and so, it contributes 2>kT/2. This principle is called the the energy equipartition 
theorem. In the sequel, we will see that it is quite intimately related to rate-distortion theory 
for quadratic distortion measures. 
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Below is a direct derivation of the equipartition theorem: 



-aX' 



r dxe-/5"-V2) 

J — oo 



num. & den. have closed forms, but we use another way: 



dp 



-A In 

dp 



^ 1 
Idln^ 



-00 
00 



1 
1 

1 _ A;T 
2 d^ '2p~^ 



d(V^x)e-"(^^ 



dwe"""'/' 



The integral is now a constant, independent of /3. 



This simple trick, that bypasses the need to calculate integrals, can easily be extended in 
two directions at least (exercise): 



Let a; e IR" and let Si Ax, where A is a n x n positive definite matrix. This 

corresponds to a physical system with a quadratic Hamiltonian, which includes also 
interactions between pairs (e.g.. Harmonic oscillators or springs, which are coupled 
because they are tied to one another). It turns out that here, regardless of A, we get: 



{E{X)) = {\x-AX^^n.^^. 



(56) 



Back to the case of a scalar x, but suppose now a more general power-law Hamiltoinan, 
S{x) — a\x\^ . In this case, we get 

kT 



(Six)) = {a\Xf) 



e 



(57) 



Moreover, if Yimy.^±^xe ^^^^^ — for all (3 > 0, and we denote S'{x) = dS{x)/dx, 
then 

{X ■ S'{X)) = kT. (58) 

It is easy to see that the earlier power-law result is obtained as a special case of this, 
as E'{x) — a9\x\^~^sgn{x) in this case. 
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Example/Exercise - Ideal gas with gravitation: Let 

E{x) = ^2m^^' + ^^^^ 

The average kinetic energy of each particle is 3kT/2, as said before. The contribution of the 
average potential energy is kT (one degree of freedom with ^ = 1). Thus, the total is bkT/2, 
where 60% come from kinetic energy and 40% come from potential energy, universally, that 
is, independent of T, m, and g. □ 

2.6 The Grand— Canonical Ensemble (Optional) 

Looking a bit back, then a brief summary of what we have done thus far, is the following: 
we started off with the microcanonical ensemble, which was very rcstricitve in the sense 
that the energy was held strictly fixed to the value of the number of particles was held 
strictly fixed to the value of n, and at least in the example of a gas, the volume was also held 
strictly fixed to a certain value V . In the passage from the microcanonical ensemble to the 
canonical one, we slightly relaxed the first of these parameters - E: Rather than insisting on 
a fixed value of E, we allowed energy to be exchanged back and forth with the environment, 
and thereby to slightly fiuctuate (for large n) around a certain average value, which was 
controlled by temperature, or equivalently, by the choice of ^. This was done while keeping 
in mind that the total energy of both system and heat bath must be kept fixed, by the 
law of energy conservation, which allowed us to look at the combined system as an isolated 
one, thus obeying the microcanonical ensemble. We then had a one-to-one correspondence 
between the extensive quantity E and the intensive variable /3, that adjusted its average 
value. But the other extensive variables, like n and V were still kept strictly fixed. 

It turns out, that wc can continue in this spirit, and 'relax' also either one of the other 
variables n V (but not both at the same time), allowing it to fiuctuate around a typical 
average value, and controlling it by a corresponding intensive variable. Like E, both n and 
V are also subjected to conservation laws when the combined system is considered. Each 
one of these relaxations, leads to a new ensemble in addition to the microcanonical and 
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the canonical ensembles that we have already seen. In the case where it is the variable n 
that is allowed to be flexible, this ensemble is called the grand-canonical ensemble. In the 
case where it is the variable V ^ this is called the Gibbs ensemble. And there are, of course, 
additional ensembles based on this principle, depending on what kind of the physical sytem 
is under discussion. We will not delve into all of them here because this not a course in 
physics, after all. We will describe, however, in some level of detail the grand-canonical 
ensemble. 

The fundamental idea is essentially the very same as the one we used to derive the 
canonical ensemble, we just extend it a little bit: Let us get back to our (relatively small) 
subsystem, which is in contact with a heat bath, and this time, let us allow this subsystem 
to exchange with the heat bath, not only energy, but also matter, i.e., particles. The heat 
bath consists of a huge reservoir of energy and particles. The total energy is Eq and the 
total number of particles is Uq. Suppose that we can calculate the density of states of the 
heat bath as function of both its energy E' and amount of particles n', call it Qhb{E', n'). A 
microstate now is a combnination (cc, n), where n is the (variable) number of particles in our 
subsystem and x is as before for a given n. Prom the same considerations as before, whenever 
our subsystem is in state (a;, n), the heat bath can be in any one of VLhe^Eq — £{x),no — n) 
microstates of its own. Thus, owing to the microcanonical ensemble. 



P{x,n) oc 



Q.hb{Eo -£{x),no- n) 



exp{SHB{Eo -S{x),no - n)/k} 




(60) 



where we have now defined the chemical potential fi (of the heat bath) as: 




(61) 



Thus, we now have the grand-canonical distribution: 

^/3[lin-£{X)] 



(62) 
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where the denominator is called the grand partition function: 

oo oo 

/.) ^ ^ e^'^- ^'''^""^ = E ^''""ZM- (63) 

n=0 X n=0 

It is sometimes convenient to change variables and to define z = e^'^ (which is called the 
fugacity) and then, define 

oo 

= (64) 

n=0 

This notation emphasizes the fact that for a given /3, H(z) is actually the z-transform of 
the sequence Z^. A natural way to think about P{x,n) is as P{n) ■ P{x\n), where P{n) is 
proportional to z^Zn{^) and P{x\n) corresponds to the canonical ensemble as before. 

Using the grand partition function, it is now easy to obtain moments of the RV n. For 

example, the first moment is: 

,._ T.nnz-Z^{P) _ ainS(/3,.) 

^^'t::^^~ ■ ^^^^ 

Thus, we have replaced the fixed number of particles n by a random number of particles, 
which concentrates around an average controlled by the parameter or equivalently, z. 
The dominant value of n is the one that maximizes the product z^Zn{f3), or equivalently, 
P/in + In Zn{P). Thus, In 5 is related to In Z„ by another kind of a Legendre transform. 

When two systems, with total energy Eq and a total number of particles no, are brought 
into contact, allowing both energy and matter exchange, then the dominant combined states 
are those for which Jli(£'i, rii) • 0:2{Eq — E-i, uq — rii), or equivalently, Si{Ei,ni) + S2{Eo — 
Ei,no — ni), is maximum. By equating to zero the partial derivatives w.r.t. both Ei and ni, 
we find that in equilibrium both the temperatures Ti and T2 are the same and the chemical 
potentials /ii and 112 are the same. 

Finally, I would hke to point out that beyond the obvious physical significance of the 

grand-canonical ensemble, sometimes it proves useful to work with it from the reason of 
pure mathematical convenience. This is shown in the following example. 

Example - Quantum Statistics. Consider an ensemble of indistinguishable particles, each 
one of which may be in a certain quantum state labeled by 1, 2, . . . , r, . . .. Associated with 
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quantum state number r, there is an energy e^. Thus, if there are particles in each state 
r, the total energy is '^j.nr€r, and so, the canonical partition function is: 

ZM= exp{-/3^n,e,}. (66) 

The constraint rir = n, which accounts for the fact that the total number of particles 
must be n, causes an extremely severe headache in the calculation. However, if we pass to 
the grand-canonical ensemble, things becomes extremely easy: 

n>0 Tl: J]^n^=n r 

^ 5^ 5^ • • • ^^"^"^ exp{-/3 Y ^^rer} 

ni>0n2>0 r 

^ I] • • • n exp{-/3n,e,} 

ni>0n2>0 r>l 
r>l nr>0 

In the case where rir is unlimited {Bose-Einstein particles, or Bosons), each factor indexed 
by r is clearly a geometric series, resulting in S = Hr [-'-/(-'- ~ ze~^^'')]. In the case where 
no quantum state can be populated by more than one particle, owing to Pauli's exclusion 
principle {Fermi-Dirac particles, or Fermions) , each factor in the product contains two terms 
only, pertaining to = 0, 1, and the result is 5 = Ylri^ + ze~^'^''). In both cases, this is 
fairly simple. Having computed E{f3,z), we can in principle, return to Zn{f3) by applying 
the inverse 2;-transform. We will get back to this in the sequel. 

2.7 Gibbs' Inequality, the 2nd Law, and the Data Processing Thm 

While the laws of physics draw the boundaries between the possible and the impossible in 
Nature, the coding theorems of information theory, or more precisely, their converses, draw 
the boundaries between the possible and the impossible in coded communication systems 
and data processing. Are there any relationships between these two facts? 

We are now going to demonstrate that there are some indications that the answer to 
this question is affirmative. In particular, we are going to see that there is an intimate 
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relationship between the second law of thermodynamics and the data processing theorem 
(DPT), asserting that if X ^ U ^ V is a Markov chain, then I{X; U) > I{X; V). The 
reason for focusing our attention on the DPT is that it is actually the most fundamental 
inequality that supports most (if not all) proofs of converse theorems in IT. Here are just a 
few points that make this quite clear. 

1. Lossy/lossless source coding: Consider a source vector — {Ui, . . .Un) compressed 
into a bitstream X" — {Xi, . . . , X^) from which the decoder generates a reproduction 

= (Vi,...,Vn) with distortion E{d{Ui,Vi}} < ND. Then, by the DPT, 

I{U^;V^) < /(X";X") = if(X"), where I{U^;V^) is further lower bounded by 
NR{D) and H{X"') < n, which together lead to the converse to the lossy data com- 
pression theorem, asserting that the compression ratio n/N cannot be less than R{D). 
The case of lossless compression is obtained as a special case where D — 0. 

2. Channel coding under bit error probability: Let — {Ui, . . . C/jv) be drawn from the 
binary symmetric course (BSS), designating M — 2^ equiprobable messages of length 
N. The encoder maps C/^ into a channel input vector X", which in turn, is sent across 
the channel. The receiver observes a noisy version of X", and decodes the message 
as . Let Pb = jj: ^f^i Pr{Vi 7^ Ui} designate the bit error probability. Then, 
by the DPT, I{U^;V^) < /(X";y"), where 7(X";y") is further upper bounded 
by nC, C being the channel capacity, and I(U^;V^) = H(U^) - H(U^\V^) > 
N - Y.tiH{Um > N - E>2(Pr{V^^ ^ Ui}) > N[l - h2{Pb)]. Thus, for to 
vanish, the coding rate, N/n should not exceed C. 

3. Channel coding under block error probability - Fano 's inequality: Same as in the pre- 
vious item, except that the error performance is the block error probability Pb = 
Pt{V^ ^ U^}. This, time H{U^\V^), which is identical to H{U^,E\V^), with 
E = I{V^ ^ U^}, is decomposed as H{E\V^) + H{U^\V^, E), where the first term 
is upper bounded by 1 and the second term is upper bounded by Pb log(2-^ — 1) < NPb, 
owing to the fact that the maximum of H{U^\V^, E — 1) is obtained when U^ is dis- 
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tributed uniformly over all V'^ ^ . Putting these facts all together, we obtain 
Fano's inequality Pb > 1 — 1/n — C/R, where R — N/n is the coding rate. Thus, the 
DPT directly supports Fano's inequality, which in turn is the main tool for proving 
converses to channel coding theorems in a large variety of communication situations, 
including network configurations. 

4. Joint source-channel coding and the separation principle: In a joint source-channel 
situation, where the source vector is mapped to a channel input vector X" and 
the channel output vector is decoded into a reconstruction , the DPT gives 
rise to the chain of inequalities NR{D) < I{U^; V^) < /(X"; Y"") < nC, which is the 
converse to the joint source-channel coding theorem, whose direct part can be achieved 
by separate source- and channel coding. Items 1 and 2 above are special cases of this. 

5. Conditioning reduces entropy: Perhaps even more often than the term "data process- 
ing theorem" can be found as part of a proof of a converse theorem, one encounters 
an equivalent of this theorem under the slogan "conditioning reduces entropy". This 
in turn is part of virtually every converse proof in the literature. Indeed, if {X, U, V) 
is a triple of RV's, then this statement means that H{X\V) > H{X\U,V). If, in 
addition, X — )> f/ V is a Markov chain, then H{X\U,V) = H{X\U), and so, 
H[X\V) > H{X\U), which in turn is equivalent to the more customary form of the 
DPT, I{X;U) > I{X;V), obtained by subtracting H{X) from both sides of the en- 
tropy inequality. In fact, as we shall see shortly, it is this entropy inequality that 
lends itself more naturally to a physical interpretation. Moreover, we can think of 
the conditioning-reduces-entropy inequality as another form of the DPT even in the 
absence of the aforementioned Markov condition, because X ^ {U,V) ^ V is always 
a Markov chain. 

Turning now to the physics point of view, consider a system which may have two possibile 
Hamiltonians - So{x) and Si{x). Let Zi{P), denote the partition function pertaining to 
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that is 

zm = Y.^-''"^''\ ^ = 0,1. 

X 

The Gibbs' inequality asserts that 

lnZi(/3) > \nZo{^3) + (5{So{X)-S^{X))o 



(68) 



(69) 



where {■)o denotes averaging w.r.t. Pq ^ the canonical distribution pertaining the Hamihonian 
£o{-)- Equivalently, this inequahty can be presented as follows: 



{£^iX)-£oiX))o> 



\ ln^i(/3)l 




r ln^o(/3)l 


[ /3 J 




[ /3 J 



Fi — Fq, 



(70) 



where Fi is the free energy pertaining to the canonical ensemble of Si, i = 0,1. 

This inequality is easily proved by defining an Hamiltoinan £x{x) = (1 — X)So{x) + 
\£i{x) = £o{x) + \ [Si{x) — £q{x)] and using the convexity of the corresponding log-partition 
function w.r.t. A. Specifically, let us define the partition function: 



(71) 



X 



Now, since £\{x) is afiine in A, then it is easy to show that \n Zx/dX^ > (just like this 
was done with d^ In Z(/3)/d/3^ > before) and so \nZ\{l3) is convex in A for fixed /3. It 
follows then that the curve of lnZx{P), as a function of A, must lie above the straight line 
that is tangent to this curve at A = (see Fig. [2]), that is, the graph corresponding to the 



affine function In Zo(/3) + A 



A=0 



. In particular, setting A = 1, we get: 



InZi(A) >lnZo(/3) + 



d\nZx{l3) 



and the second term is: 



d\nZx{l3) 



Thus, we have obtained 



dX 



(72) 



A=0 



(3 {SoiX) - £,iX)) 



' 



(73) 



In 



> In 



X 



-fSSoiX) 



+ (3{£oiX)-£,{X)),. 



(74) 
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straight line - tangent at 
A = 0. 



1 



-A 



Figure 2: The function In Zx(l3) is convex in A and hence hes above its tangent at the origin. 

and the proof is complete. In fact, the l.h.s. minus the r.h.s. is nothing but D{Pq\\Pi), where 
Pi is the B-G distribution pertaining to i = 0,1. 

We now offer a possible physical interpretation to the Gibbs' inequality: Imagine that a 

system with Hamiltoinan Sq{x) is in equilibrium for all t < 0, but then, at time t = 0, the 

Hamitonian changes abruptly from the Sq{x) to Si{x) (e.g., by suddenly applying a force on 

the system), which means that if the system is found at state x at time t = 0, additional 

energy of = Si{x) — So{x) is suddenly 'injected' into the system. This additional energy 

can be thought of as work performed on the system, or as supplementary potential energy. 

Since this passage between £q and £i is abrupt, the average of W should be taken w.r.t. Pq, 

as the state x does not change instantaneously. This average is exactly what we have at the 

left-hand side eq. (*). The Gibbs inequality tells us then that this average work is at least 

as large as AF = Fi — Fq, the increase in free energy^ The difference {W)q — AF is due 

to the irreversible nature of the abrupt energy injection, and this irreversibility means an 

increase of the total entropy of the system and its environment, and so, the Gibbs' inequality 

^This is related to the interpretation of the free-energy difference AF ^ Fi — Fq as being the maximum 
amount of work in an isothermal process. 
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is, in fact, a version of the second law of thermodynamics|l| This excess work beyond the 
free-energy increase, {W)q — AF, which can be thought of as the "dissipated work," can 
easily shown (exercise) to be equal to kT ■ D (Poll -Pi), where Pq and Pi are the canonical 
distributions pertaining to Sq and £i, respectively. Thus, the divergence is given yet another 
physical significance. 

Now, let us see how the Gibbs' inequality is related to the DPT. Consider a triple of 
random variables {X, U, V) which form a Markov chain X ^ U ^ V . The DPT asserts 
that I{X; U) > I{X; V). We can obtain the DPT as a special case of the Gibbs inequal- 
ity as follows: For a given realization {u,v) of the random variables {U,V), consider the 
Hamiltonians 

£q(x) = —In P{x\u) = —hi P{x\u,v) (75) 

and 

£i{x) = -\nP{x\v). (76) 
Let us also set (3 = 1. Thus, for a given {u, v): 

{W)o = {Si{X)-So{X))o = J^P(£c|w,'y)[lnP(a;|w)-lnP(a;|v)] = H{X\V = v)-H{X\U = 

X 

(77) 

and after further averaging w.r.t. {U, V), the average work becomes H{X\V) — H{X\U) = 
I{X; U) — I{X; V). Concerning the free energies, we have 

Zo{(3 = 1) = ^exp{-l ■ [-\nP{x\u,v)]} = ^ P(a;|w, -u) = 1 (78) 

X X 

and similarly, 

Zi(/3 = l) = 5^P(a;|t;) = l (79) 

X 

^ From a more general physical perspective, the Jarzynski equality tells that under certain conditions 
on the test system and the heat bath, and given any protocol {X{t)} of changing the control variable A (of 
£\{x)), the work W applied to the system is a RV which satisfies (e~^^) = e~^^^ . By Jensen's inequality, 
^g-/3VK^ is lower bounded by e~^^^\ and so, we obtain [W) > AF (which is known as the minimum work 
principle), now in more generality than in the Gibbs' inequality, which is limited to the case where X{t) is a 
step function. At the other extreme, when X(t) changes very slowly, corresponding to a reversible process, 
W approaches determinism, and then Jensen's inequality becomes tight, which then gives (in the limit) 
W = AF with no increase in entropy. 
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which means that Fq — Fi — 0, and so AF = as well. So by the Gibbs inequahty, the 
average work /(X; U) — I{X; V) cannot be smaller than the free-energy difference, which 
in this case vanishes, namely, I{X; U) — I{X; V) > 0, which is the DPT. Note that in 
this case, there is a maximum degree of irreversibility: The identity I{X] U) — I{X] V) — 
H{X\V) - H{X\U) means that whole work W = I{X]U) - I{X]V) goes for entropy 
increase SiT — SqT — H{X\V) ■ 1 — H{X\U) • 1, whereas the free energy remains unchanged, 
as mentioned earher. Note that the Jarzynski formula (cf. last footnote) holds in this special 
case, i.e., (e-^-^) = e"^'^^' = 1. 

The difference between I{X; U) and I{X] V), which accounts for the rate loss in any 
suboptimal coded communication system, is then given the meaning of irreversibility and 
entropy production in the corresponding physical system. Optimum (or nearly optimum) 
communication systems arc corresponding to quasistatic isothermal processes, where the full 
free energy is exploited and no work is dissipated (or no work is carried out at all, in the first 
place). In other words, had there been a communication system that violated the converse to 
the source/channel coding theorem, one could have created a corresponding physical system 
that violates the second law of thermodynamics, and this, of course, cannot be true. 

2.8 Large Deviations Theory and Physics of Information Measures 

As I said in the Intro, large deviations theory, the branch of probabihty theory that deals 
with exponential decay rates of probabilities of rare events, has strong relations to IT, which 
we have already seen in the IT course through the eye glasses of the method of types and 
Sanov's theorem. On the other hand, large deviations theory has also a strong connection 
to statistical mechanics, as we are going to see shortly. Therefore, one of the links between 
IT and statistical mechanics goes through rate functions of large deviations theory, or more 
concretely, ChernofT bounds. This topic is based on the paper: N. Merhav, "An identity of 
Chernoff bounds with an interpretation in statistical physics and applications in information 
theory," IEEE Trans. Inform. Theory, vol. 54, no. 8, pp. 3710-3721, August 2008. 

Let us begin with a very simple question: We have a bunch of i.i.d. RV's Xi, X2, . . . and 
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a certain real function £{x). How fast does the probability of the event 

n 

Y,S{Xi)<nE^ 



i=l 



decay as n grows without bound, assuming that Eq < {£{X)) (so that this would be a rare 
event)? One way to handle this problem, at least in the finite alphabet case, is the method 
of types. Another method is the Chernoff bound: 

Pr < ^ £{Xi) < uEq > = EX < ^ £{Xi) < hEq > X(-) denoting the indicator function 



i=l 



nEo-^S{Xi 



i=l 



^ V /3 > : I{Z <a} <e 



P{a-Z) 



i=l 



= e'^"^°i;|ncxp{-^£:(x,)} 

= e'^"^°[i;exp{-/3£(Xi)}]" 

= exp{n[^£;o + ln£;exp{-^£:(Xi)}]} 

As this bound applies for every /5 > 0, the tightest bound of this family is obtained by 
minimizing the r.h.s. over ^, which yields the exponential rate function: 



where 



and 



E(£;o) = min[/3£;o + (/)(^)], 

/3>0 



0(^) = lnZ(^) 



(80) 



(81) 



(82) 



Rings a bell? Note that Z[j3) here differs from the partition function that we have encoun- 
tered thus far only slighlty: the Boltzmann exponentials are weighed by {p{x)} which are 
independent of ^. But this is not a crucial difference: one can imagine a physical system 
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where each microstate x is actually a representative of a bunch of more refined microstates 
{x'}, whose number is proportional to p{x) and which all have the same energy as x, that 
is, £{x') = S{x). In the domain of the more refined system, Z{f3) is (up to a constant) a 
non-weighted sum of exponentials, as it should be. More precisely, if p{x) is (or can be 
approximated by) a rational number N{x)/N, where N is independent of x, then imagine 
that each x gives rise to N{x) microstates {x'} with the same energy as x, so that 



and we are back to an ordinary, non-weighted partition function, upto the constant 1/N, 
which is absolutely immaterial. 

To summarize what we have seen thus far: the exponential rate function is given by the 
Lcgcndrc transform of the log-moment generating function. The Chernoff parameter /3 to 
be optimized plays the role of the equilibrium temperature pertaining to energy Eq. 

Consider next what happens when p{x) is itself a B-G distribution with Hamiltonian 

S{x) at a certain inverse temperature that is 




(83) 



X 



x' 



e 



■0ie{x) 



p{x) 



(84) 



with 




(85) 



X 



In this case, we have 




(86) 



Thus, 



min[^£;o + lnC(A + ^)]-lnC(A) 

mm[{l3 + Pi)Eo + In ((A + P)] - lnC(/3i) - A^o 

^m[l3Eo + lnC(/3)] - lnC(/3i) - /^.E^ 

mm[/3Eo + In C(^)] - [In C(A) + PiEi] + /3i{Ei - Eo) 
p>pi 
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where Ei is the energy corresponding to i.e., Ei is such that 



a{E,) ^ 



min[/3Ei + lnC(/3)] 



(87) 



is achieved by /3 = Thus, the second bracketted term of the right-most side of the last 
chain is exactly a{Ei), as defined. If we now assume that Eq < Ei, which is reasonable, 
because Ei is the average of S{X) under and we are assuming that we are dealing with 
a rare event where Eq < {S{X)). In this case, the achiever /3o of cr(i?o) must be larger than 
Pi anyway, and so, the first bracketted term on the right-most side of the last chain agrees 
with a{Eo). We have obtained then that the exponential decay rate (the rate function) is 
given by 



Note that / > thanks to the fact that cr(-) is concave. It has a simple graphical intepretation 
as the height difference, as seen at the point E = Eq, between the tangent to the curve cr{E) 
at E = El and the function a{E) itself (see Fig. [3]). 



S(Eo) = a{Ei) - a{Eo) - /3i{Ei - E^). 



(88) 



a{E) 




E 



Figure 3: Graphical interpretation of the LD rate function I. 
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Another look is the following: 

= /3i(Fo-Fi) 
= D{Pf,,\\Pp,) 

= min{D(g||P^J : Eq£{X) < Eq} ^ exercise 

The last line is exactly what we would have obtained using the method of types. This means 
that the dominant instance of the large deviations event under discussion pertains to thermal 
equilibrium (minimum free energy) complying with the constraint (s) dictated by this event. 
This will also be the motive of the forthcoming results. 

Exercise: What happens if p{x) is B-G with an Hamiltonian £(■), different from the one 
of the LD event? □ 

Let us now see how this discussion relates to very fundamental information measures, like 
the rate-distortion function and channel capacity. To this end, let us first slightly extend 
the above Chernoff bound. Assume that in addition to the RV's Xi, . . . , there is also a 
deterministic sequence of the same length, yi, . . . ,yn, where each yi takes on values in a finite 
alphabet y. Suppose also that the asymptotic regime is such that as n grows without bound, 
the relative frequencies Y17=i ^iVi ~ y}}yey converge to certain probabilities {q{y)}yey- 
Furthermore, the Xj's are still independent, but they are no longer necessarily identically 
distributed: each one of them is governed by p{xi\yi), that is, p{x\y) = YYi=iPi^i\yi)- Now, 
the question is how does the exponential rate function behave if we look at the event 

n 

Y,^{Xi,yi)<nEo (89) 

where £{x,y) is a given 'Hamiltonian'. What is the motivation for this question? Where 
and when do we encounter such a problem? 

Well, there are many examples (cf. the above mentioned paper), but here are two very 
classical ones, where rate functions of LD events are directly related to very important 
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information measures. In both examples, the distributions p{-\y) are actually the same for 
all y & y (namely, {Xi} are again i.i.d.). 

• Rate-distortion coding. Consider the good old problem of lossy compression with a 
randomly selected code. Let y = {yi, . . . ,yn) be a given source sequence, typical to 
Q = {Q{y)y y y} (non-typical sequences are not important). Now, let us randomly 
select e"^ codebook vectors {X{i)} according to p{x) — Y[i=iP{^i)- Here is how the 
direct part of the source coding theorem essentially works: We first ask ourselves what 
is the probability that a single randomly selected codeword X — (Xi, . . . , X„) would 
happen to fall at distance < nD from y, i.e., what is the exponential rate of the 
probability of the event 

n 

J2d{Xi,yi)<nD7 (90) 

i=l 

The answer is that it is exponentially about e~"^(^), and that's why we need slightly 
more than one over this number, namely, e"*""^*^^^ times to repeat this 'experiment' in 
order to see at least one 'success', which means being able to encode y within distortion 
D. So this is clearly an instance of the above problem, where £ — d and Eq — D. 

• Channel coding. In complete duality, consider the classical channel coding problem, for 
a discrete memoryless channel (DMC), using a randomly selected code. Again, we have 

a code of size e"'^, where each codeword is chosen independently according to 'p{x) = 
YYi=iP{xi). Let y the channel output vector, which is (with very high probabaility) , 
typical to Q — {q{y), y e y}, where q{y) — '^j.p{x)W{y\x), W being the single-letter 
transition probability matrix of the DMC. Consider a (capacity-achieving) threshold 
decoder which selects the unique codeword that obeys 

n 

J2[-\nW(yi\Xi)]<n[H{Y\X) + e] e>0 (91) 

i=l 

and declares an error whenever no such codeword exists or when there is more than 
one such codeword. Now, in the classical proof of the direct part of the channel coding 
problem, we first ask ourselves: what is the probability that an independently selected 
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codeword (and hence not the one transmitted) X will pass this threshold? The answer 
turns out to be exponentially e~"^, and hence we can randomly select up to slightly less 

than one over this number, namely, e"'""''" codewords, before wc start to see incorrect 
codewords that pass the threshold. Again, this is clearly an instance of our problem 
with S{x, y) = - In W{y\x) and Eo = H{Y\X) + e. 

Equipped with these two motivating examples, let us get back to the generic problem we 
formalized, and see what happens. Once this has been done, we shall return to the examples. 
There are (at least) two different ways to address the problem using Chernoff bounds, and 
they lead to two seemingly different expressions, but since the Chernoff bounding technique 
gives the correct exponential behavior, these two expressions must agree. This identity 
between the two expressions will have a physical intepretation, as we shall see. 

The first approach is a direct extension of what we did before: 




Ey = expectation under p{-\y) 



= e^"^° n ^M-mX, y)}f'^ n{y) ^ num. of {y, = y} 



\ L yey xex 

and so, the resulting rate function is given by 





(92) 



where 




(93) 
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In the rate-distortion example, this tells us that 



R(D) — — min 

/3>0 



+ J]?(y)lnJ]p(x)e-'^''(^'^) 



(94) 



This is a well-known parametric representation of R{D), which can be obtained via a different 
route (see, e.g.. Gray's book Source Coding Theory), where the minimizing /3 is known to 
have the graphical interpretation of the local negative slope (or derivative) of the curve of 
R{D). In the case of channel capacity, we obtain in a similar manner: 



C 



— mm 

/9>0 



— mm 

/9>0 



I3H{Y\X) + J2 <l{y) In ^p(a;)e-'5[-^'^^(^l 



x)] 



y&y 



/3H{Y\X) + J2 liy) In J]p(x)iy'^(y|x) 
yey xeX 



Exercise: Show that for channel capacity, the minimizing (5 is always (5* = 1. □ 
The other route is to handle each y & y separately: First, observe that 

n 

^£(x„yo = E E ^(^-^)' 

1=1 yey i: yi=y 



(95) 



where now, in each partial sum over {i : yi — y}, we have i.i.d. RV's. The event 
Yl^=i ^{^h Vi) < "^Eq can then be thought of as the union of all intersections 



flj E S{X,,y)<n{y)EA 
yey Ki- yi=y ) 



(96) 



where the union is across all "possible partial energy allocations" {Ey} which satisfy <l{y)Ey 
Eq. Note that at least when the Xj's take values on a finite alphabet, each partial sum 
X^i- y^=y^{^^^y) can take only a polynomial number of values in n{y) (why?), and so, it is 
sufficient to 'sample' the space of {Ey} by polynomially many vectors in order to cover all 



< 
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possible instances of the event under discussion (see more details in the paper). Thus, 

{n 
J2S{Xi,yi)<nEo 



q{y)Ey<Eo} yey Ki: yi=y 



= Pr 

{Ey: Y,yi{y)Ey<Eo}yey U 



yey yi=y 



— max TT exp < n(y) mini By Ey + In Zy 

{Ey. Eyi(y)Ey<Eo}^^^ ^\ ^^^,>oL^^ ^ ^ 

= exp < n ■ max > q{y)TiJEy) 

I {Ey. T.y^{y)Ey<E,}j^^^^> 



where we have defined 



^y{Ey) = min [PyEy + In . (97) 

PydV 



We therefore arrived at an alternative expression of the rate function, which is 

v^f^P^<riE^(^)^^(^^)- (98) 
{Ey- T,yg{y)Ey<EQ}^ 

Since the two expressions must agree, we got the following identity: 



E(Eo) = max{E^, j:,q(y)Ey<Eo} Y^y^y Qiy)^yiEy) 



A few comments: 

1. In the paper there is also a direct proof of this identity, without relying on Chernoff bound 
considerations. 

2. This identity accounts for a certain generalized concavity property of the entropy function. 
Had all the Ej^(-)'s been the same function, then this would have been the ordinary concavity 
property. What makes it interesting is that it continues to hold for different Ej,(-)'s too. 
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3. The l.h.s. of this identity is defined by minimization over one parameter only - the inverse 
temperature /3. On the other hand, on the r.h.s. we have a separate inverse temperature 
for every because each T,y{-) is defined as a separate minimization problem with its own 
Py. Stated differently, the l.h.s. is the minimum of a sum, whereas in the r.h.s., for given 
{Ey}, we have the sum of minima. When do these two things agree? The answer is that 
it happens if all minimizers {(3*} happen to be the same. But (3* depends on Ey. So what 
happens is that the {Ey} (of the outer maximization problem) are such that the /3* would 
all be the same, and would agree also with the /3* of S(i?o)- To see why this is true, consider 
the following chain of inequalities: 

max q(y)Tjy(Ey) 

y y 



{Ey: T.y'}{y)Ey<Eo}^ Py 



y 



< 



max V q(y) \f3*Ey + In ZJ(3*)] where /3* achieves S(Eo) 



{Ey: T.yQ{y)Ey<Eo} 



y 



< 



_ max [/3*Eo + q{y) In Zy{(3*)] because q{y)Ey < Eq 

{Ey: T.yl{y)Ey<EQ} ^ 

y 

= f3*Eo + '^q{y)lnZy{P*) the bracketted expression no longer depends on {Ey} 

y 

= S(Eo). 

Both inequalities become equalities if {Ey} would be allocated such thatj^ (i) Yyqiy)Ey = 
Eq and (ii) l3*{Ey) = 13* for all y. Since the /3's have the meaning of inverse temperatures, 
what we have here is thermal equilibrium: Consider a bunch of |3^| subsystems, each one 
of n{y) particles and Hamiltonian £{x, y) indexed by y. If all these subsystems are thermally 
separated, each one with energy Ey, then the total entropy per particle is 'Yy(l{y)'^y{Ey). 
The above identity tells us then what happens when all these systems are brought into 
thermal contact with one another: The total energy per particle Eq is split among the 
different subsystems in a way that all temperatures become the same - thermal equilibrium. 
It follows then that the dominant instance of the LD event is the one where the contributions 



^Exercise: show that there exists an energy allocation {Ey} that satisfies both (i) and (ii) at the same 
time. 
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of each y, to the partial sum of energies, would correspond to equilibrium. In the rate- 
distortion example, this characterizes how much distortion each source symbol contributes 
typically. 

Now, let us look a bit more closely on the rate-distortion function: 



R(D) — — min 

/3>0 



(99) 



As said, the Chernoff parameter (3 has the meaning of inverse temperature. The inverse 
temperature P required to 'tune' the expected distortion (internal energy) to D, is the 
solution to the equation 

^--^H ^(^) In J]p(a;)e-^''(^'^) (100) 
y X 

or equivalently. 

The Legendre transform relation between the log-partition function and R{D) induces a one- 
to-one mapping between D and (3 which is defined by the above equation. To emphasize 
this dependency, we henceforth denote the value of corresponding to a given /3, by D^. 
This expected distortion is defined w.r.t. the probability distribution: 

P,{x,y) = g(y) • PMv) - Q{y) ■ ,yy (102) 

On substituting instead of D in the expression of R{D), we have 

-i?(D^) = /3D^ + ^g(|/)ln^p(a:)e-'^'^(-'^). (103) 

y X 

Note that R{Dp) can be represented in an integral form as follows: 



P ■ dD^, (104) 

Do 
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where Dq — ^xyPi^)^iy)^i^^y) value of D corresponsing to ^ = 0, and for which 

Rq{D) = 0, This is exactly analogous to the thermodynamic equation S = J dQ/T (following 
from 1/T = d5'/dQ), that builds up the entropy from the cumulative heat. Note that the last 
equation, in its differential form, reads dR{D^) = —(3dDj3, or (3 = —R'{Dp), which means 
that P is indeed the negative local slope of the rate-distortion curve R{D). Returning to 
the integration variable ^, we have: 

dD, 



R{Dp) = - / d/3-/3 
Jo 



d(5 

= Y.^{y) j d^-^-Var^{d(X,y)|y = y} 
y 

= / dp- ^■m.mseAd{X,Y)\Y} 
Jo 

where Var^{-} and mmse^{-|-} are taken w.r.t. Pp{x,y). We have therefore introduced an 
integral representation for R{D) based on the MMSE in estimating the distortion variable 
d{X, Y) based on Y. In those cases where an exact expression for R{D) is hard to obtain, 
this opens the door to upper and lower bounds on R{D), which are based on upper and lower 
bounds on the MMSE, offered by the plethora of bounds available in estimation theory. 
Exercise: Show that Dp — Dq — Jq d^ ■ mmse p{d{X,Y)\Y}. 

Finally, a word about the high-resolution regime. The partition function of each y is 

Zy{/3) = J2p{x)e~^'^^"''^\ (105) 

X 

or, in the continuous case, 

Zy{/3) = [ dxp{x)e-^'^^''^yl (106) 

Consider the distortion measure d{x,y) — \x — y\^, where ^ > and consider a uniform 
random coding distribution over the interval [— 74,^4], supposing that it is the optimal (or 
close to optimal) one. Suppose further that we wish to work at a very small distortion level 
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D (high res), which means a large value of P (why?). Then, 



ZyW) 



1 

2A 
1 

2A 
1 

2A 



A 

+ CX3 

■oo 

+ 00 



(large /3) 

dare"^'^'" (the integral is independent of y) 



Thus, returning to the expression of R{D), let us minimize over (3 by writing the zero- 
derivative equation, which yields: 



d_ 



In 



1 

2A 



+ 00 



dxe-'^l-l' 



(107) 



but this is exactly the calculation of the (generalized) cquipartition theorem, which gives 
l/{pe) = kT/e. Now, we already said that /3 = -R'{D), and so, = -D'{R). It follows 
then that the function D{R), at this high res. limit, obeys a simple differential equation: 

D'{R) 



D{R)^ 



e 



(108) 



whose solution is 

D{R) = Doe"^^. (109) 

In the case where 6 = 2 (squared error distortion), we get that D[R) is proportional to e~^^, 
which is a well-known result in high res. quantization theory. For the Gaussian source, this 
is true for all R. 
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3 Analysis Tools and Asymptotic Methods 

3.1 Introduction 



So far we have dealt with relatively simple situations where the Hamiltonian is additive, 
the resulting B-G distribution is then i.i.d., and everything is very nice, easy, and simple. 
But this is seldom the case in reahty. Most models in physics, including those that will 
prove relevant for IT, as we shall sec in the sequel, are way more complicated, more difficult, 
but also more interesting. More often than not, they are so complicated and difficult, that 
they do not lend themselves to closed-form analysis at all. In some other cases, analysis 
is possible, but it requires some more powerful mathematical tools and techniques, which 
suggest at least some asymptotic approximations. These are tools and techniques that we 
must acquaint ourselves with. So the purpose of this part of the course is to prepare these 
tools, before we can go on to the more challenging settings that are waiting for us. 

Before diving into the technical stuff, I'll first try to give the flavor of the things I am 
going to talk about, and I believe the best way to do this is through an example. In quantum 
mechanics, as its name suggests, several physical quantites do not really take on values in 
the continuum of real numbers, but only values in a discrete set, depending on the conditions 
of the system. One such quantized physical quantity is energy (for example, the energy of 
light comes in quanta of hi/, where u is frequency). Suppose we have a system of n mobile 
particles (gas), whose energies take on discrete values, denoted eo < ei < £2 < . . .. If the 
particles were not interacting, then the partition function would have been given by 



r>Q 



ri>0r-2>0 r-„>0 L i=l ) 71: Y,r'^r=n '^'^ ^' I r 



(110) 

However, since the particles are indistinguishable, then permutations among them are not 
considered distinct physical states (see earlier discussion on the ideal gas), and so, the com- 
binatorial factor n\/ W^nr\i that counts these permutations, should be eliminated. In other 
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(113) 



words, the correct partition function should be 

ZM- Yl exp|-^^n,eJ. (Ill) 

The problem is that this partition function is hard to calculate in closed form: the headache 
is caused mostly because of the constraint = n. However, if we define a corresponding 

generating function 

?i>0 

which is like the ^-transform of {Z„(/3)}, this is easy to work with, because 
H(/3, .) ^ J] J] . . . .2:.- exp 5^n.e.| = n [E (^^"'")" 

ni>0n2>0 V r ) r \_ rir 

Splendid, but we still want to obtain Zn{P)... 
The idea is to apply the inverse z-transform: 

Ufi) - ^ I - ^ I m .)e-(»«)-d., (114) 

where 2; is a complex variable, j — and C is any clockwise closed path encircling the 

origin and entirely in the region of convergence. An exact calculation of integrals of this 
type might be difficult, in general, but often, we would be happy enough if at least we could 
identify how they behave in the thermodynamic limit of large n. 

Similar needs are frequently encountered in information-theoretic problems. One exam- 
ple is in universal source coding: Suppose we have a family of sources indexed by some 
parameter 9, say, Bernoulli with parameter 9 e [0, 1], i.e.. 



Pe(a;) = (1 - ^) xe{0,iy; n = #ofl's (115) 

When 9 is unknown, it is customary to construct a universal code as the Shannon code w.r.t. 
a certain mixture of these sources 

P(a;)= / d9w{9)P0{x)^ [ d9w{9)e^^^'^^ (116) 
JO Jg 
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where 

M^)=ln(l-^) + gln(^j; q=^. (117) 

So here again, we need to evaluate an integral of an exponential function of n (this time, on 
the real line), in order to assess the performance of this universal code. 

This is exactly the point where the first tool that we are going to study, namely, the 
saddle point method (a.k.a. the steepest descent method) enters into the picture: it gives us 
a way to assess how integrals of this kind scale as exponential functions of n, for large n. 
More generally, the saddle point method is a tool for evaluating the exponential order (plus 
2nd order behavior) of an integral of the form 

g{z)e^^^^'^(\z V is Si path in the complex plane. (US) 

We begin with the simpler case where the integration is over the real line (or a subset of 
the real line), whose corresponding asymptotic approximation method is called the Laplace 
method. The material here is taken mostly from de Bruijn's book, which appears in the 
bibliographical list. 

3.2 The Laplace Method 

Consider first an integral of the form: 

/•+0O 

F„= / e'^'^^^Mx, (119) 

J —DC 

where the function h{-) is independent of n. How does this integral behave exponentially 
for large n? Clearly, if it was a sum, like X^i^^'^S rather than an integral, and the number 
of terms was finite and independent of n, then the dominant term, e"™^''*% would have 
dictated the exponential behavior. This continues to be true even if the sum contains even 
infinitely many terms provided that the tail of this series decays sufficiently rapidly. Since 
the integral is, after all, a limit of sums, it is conceivable to expect, at least when h{-) is 
"sufficiently nice", that something of the same spirit would happen with F„, namely, that its 
exponential order would be, in analogy, eJ^^^K^) , In what follows, we are going to show this 



T 
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more rigorously, and as a bonus, we will also be able to say something about the second order 
behavior. In the above example of universal coding, this gives rise to redundancy analysis. 

We will make the following assumptions on h: 

1. h is real and continuous. 

2. h is maximum at x — and h{0) — (w.l.o.g). 

3. h{x) < Vx 7^ 0, and 36 > 0, c > s.t. \x\ > c implies h{x) < —b. 

4. The integral defining F„ converges for all sufficiently large n. W.l.o.g., let this suffi- 
ciently large n he n — 1, i.e., e'^^^Mx < oo. 

5. The derivative h'{x) exists at a certain neighborhood oi x = 0, and h"{0) < 0. Thus, 
h'{0) = 0. 

Prom these assumptions, it follows that for all S > 0, there is a positive number 'r){S) s.t. for 
all \x\ > 5, we have h{x) < —r){5). For 5 > c, this is obvious from assumption 3. li 5 < c, 
then the maximum of the continuous function h across the interval [5, c] is strictly negative. 
A similar argument applies to the interval [— c, —S]. Consider first the tails of the integral 
under discussion: 

[ e^'^^^Mx = [ dxe("-^)''(^)+''(^) 

J\x\>S J\^\>S 

< I dxe-("-^)''(^)+'^(^) 

J\x\>5 

/ + 0O 
e'^^^Mx exponentially fast 
-oo 

In other words, the tails' contribution is vanishingly small. It remains to examine the integral 
from —5 to +5, that is, the neighborhood of x = 0. In this neighborhood, we shall take the 
Taylor series expansion of h. Since /i(0) = /i'(0) = 0, then h{x) fa |/i"(0)x^. More precisely, 
for all e > 0, there is 5 > s.t. 

< ex^ y\x\ < 5. (120) 



h{x) - -h"{0)x^ 
2 
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Thus, this integral is sandwiched as follows: 

y"^^exp|^(/i"(0) -e)a;2}da; < ^^"^ e"'^(^Mx < exp |^(/i"(0) + e)^^} dx. (121) 

The right-most side is further upper bounded by 

exp{-{h"{0) + e)x^]dx (122) 

and since h"{0) < 0, then h"{0) + e = —{\h"{0) \ — e), and so, the latter is a Gaussian integral 
given by 

(123) 



i\h"{0)\ -e)n 

The left-most side of the earlier sandwich is further lower bounded by 

y"^'exp{-^(|/i"(0)| + 6)x2}dx 

exp \~{\h"{0)\ + e)x^} dx- f exp \~{\h"{0)\ + e)^') da; 
'-2 J 12 J 



+ 00 



\x\>S 

- 2Q{Sy/n{\h"{0)\ + e)) 



> 



\h"(0)\ + e)n 



27r 



\h"{0)\ + e)n 
2n 

{\h"{0)\ + e)n 

where the notation An ~ Bn means that lim„_^oo ^n/-Bn = 1- Since e and hence 5 can be 
made arbitrary small, we find that 

"+(5 
-5 



Finally, since the tails contribute an exponentially small term, which is negligible compared 
to the contribution of 0(1/ -^/n) order of the integral across [—6, +d], we get: 



+ 00 



27r 



e 



oo 
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Slightly more generally, if h is maximized at an arbitrary point x — Xq this is completely 
immaterial because an integral over the entire real line is invariant under translation of the 
integration variable. If, furthermore, the maximum h{xQ) is not necessarily zero, we can 
make it zero by decomposing h according to h{x) = h{xo) + [h{x) — h{xo)] and moving the 
first term as a constant factor of e"-'^^^°^ outside of the integral. The result would then be 



Of course, the same considerations continue to apply if F„ is defined over any finite or half- 
infinite interval that contains the maximizer x — 0, or more generally x — xq a,s an internal 
point. It should be noted, however, that if F„ is defined over a finite or semi-infinite interval 
and the maximum of h is obtained at an edge of this interval, then the derivative of h at that 
point does not necessarily vanish, and the Gaussian integration would not apply anymore. In 
this case, the local behavior around the maximum would be approximated by an exponential 
exp{—n\h' {0)\x} or exp{—n\h' {xo)\x} instead, which gives a somewhat different expression. 
However, the factor e"'*(^°\ which is the most important factor, would continue to appear. 
Normally, this will be the only term that will interest us, whereas the other factor, which 
provides the second order behavior will not be important for us. A further extension in the 
case where the maximizer is an internal point at which the derivative vanishes, is this: 



where g is another function that does not depend on n. This technique, of approximating 
an integral of a function, which is exponential in some large parameter n, by neglecting the 
tails and approximating it by a Gaussian integral around the maximum, is called the Laplace 
method of integration. 



We now expand the scope to integrals along paths in the complex plane, which are also 
encountered and even more often than one would expect (cf. the earlier example). As said. 




(126) 



\h"{xo)\n 



3.3 The Saddle Point Method 
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the extension of the Laplace integration technique to the complex case is called the saddle- 
point method or the steepest descent method, for reasons that will become apparent shortly. 
Specifically, we are now interested in an integral of the form 

F^= [ e^'^^'Mz or more generally F„ = /" g{z)e''''^'^Az (127) 
J-p Jv 

where z = x + jy is a, complex variable (j = -\/— T), and P is a certain path (or curve) in the 

complex plane, starting at some point A and ending at point B. We will focus first on the 

former integral, without the factor g. We will assume that V is fully contained in a region 

where h is analytic (differentiable as many times as we want). 

The first observation, in this case, is that the value of the integral depends actually only 
on A and -B, and not on the details of V: Consider any alternate path V from AtoB such 
that h has no singularities in the region surrounded by V\JV'. Then, the integral of e"^*^^'' 
over the closed path V\JV' (going from AtoB via V and returning to A via V) vanishes, 
which means that the integrals from A to B via V and via V are the same. This means 
that we actually have the freedom to select the integration path, as long as we do not go too 
far, to the other side of some singularity point, if there is any. This point will be important 
in our forthcoming considerations. 

An additional important observation has to do with yet another basic property of analytic 
functions: the maximum modulus theorem^ which basically tells that the modulus of an 
analytic function has no maxima. We will not prove here this theorem, but in a nutshell, 
the point is this: Let 

h{z) = u{z) + jv{z) = u{x, y) + jv{x, y), (128) 

where u and v are real functions. If h is analytic, the following relationships (a.k.a. the 
Cauchy-Riemann conditions)o between the partial derivatives of u and v must hold: 

du dv _ du dv (129) 
dx dy^ dy dx 



^This is related to the fact that for the derivative f'{z) to exist, it should be independent of the direction 
at which z is perturbed, whether it is, e.g., the horizontal or the vertical direction, i.e., /'(z) = lmis^a[f {z + 
S) — f{z)]/S — lini5_>o[/(z + jS) — f{z)]/{j6), where 6 goes to zero along the reals. 
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Taking the second order partial derivative of u: 

d'^u d'^v d'^v d'^u 



(130) 



dx^ dxdy dydx dy^ 

where the first equality is due to the first Cauchy-Riemann condition and the tliird equality 

is due to the second Cauchy-Riemann condition. Equivalently, 

d'^u ^ d'^u _ ^ 
dx"^ dy^ ' 

which is the Laplace equation. This means, among other things, that no point at which 
du/dx — du/dy — can be a local maximum (or a local minimum) of because if it is 
a local maximum in the x-direction, in which case, d'^ujdx^ < 0, then d^u/dy"^ must be 
positive, which makes it a local minimum in the y-direction, and vice versa. In other words, 
every point of zero partial derivatives of u must be a saddle point. This discussion applies 
now to the modulus of the integrand e"'*^^^ because 



exp{nh{z)} 



exp[nRe{/i(^)}] = e""^^). (132) 



Of course, if h'{z) = at some z = Zq, then ^'(^o) = too, and then zq is a saddle point of 
xhus, zero-derivative points of h are saddle points. 

Another way to see this is the following: Given a complex analytic function f{z), we 
argue that the average of / over a circle always agrees with its value at the center of this 
circle. Specifically, consider the circle of radius R centered at Zq, i.e., z = Zo + Re^^ . Then, 

1 r / (^0 + Re^^) jRe^^de 



Re^^ 

1 r f {zq + Re^') d {zq + Re^') 

' ^ = f{z«). (133) 



27rj Jz=zo+ReJO Z- Zq 



and so, 

1/(^0)1 <^£ 



/ (^0 + Re'') 



d9 (134) 



which means that |/(-2o)| cannot be strictly larger than all \f{z)\ in any neighborhood (an 
arbitrary radius R) of Zq. Now, apply this fact to f{z) — e"^^^\ 
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Equipped with this background, let us return to our integral F^. Since we have the 
freedom to choose the path V, suppose that we can find one which passes through a saddle 
point zq (hence the name of the method) and that max^gp |e"'^(^)| is attained at zq. We 
expect then, that similarly as in the Laplace method, the integral would be dominated by 
g"./i(2o)^ Of course, such a path would be fine only if it crosses the saddle point zq at a 
direction w.r.t. which zq is a local maximum of |e"''*^^''|, or equivalently, of u{z). Moreover, in 
order to apply our earlier results of the Laplace method, we will find it convenient to draw 
V such that any point z in the vicinity of zq, where in the Taylor expansion is: 

h{z) ^ h{zo) + ^h"{zo){z - zoY (recall that /i'(zo) = 0.) (135) 

the second term, ^h"{zQ){z — Zq)'^ is purely real and negative, and then it behaves locally 
as a negative parabola, just like in the Laplace case. This means that 

aTg{h"{zo)} + 2arg(2 - zo) = tt (136) 

or equivalently 

arg(2; - zq) = = d. (137) 

Namely, V should cross zo in the direction 6. This direction is called the axis of zq, and 
it can be shown to be the direction of steepest descent from the peak at zq (hence the 
name) |^ 

So pictorially, what we are going to do is choose a path V from A to B, which will be 
composed of three parts (see Fig. H]): The parts A ^ A' and B' B are quite arbitrary as 
they constitute the tail of the integral. The part from A' to B', in the vicinity of zq, is a 
straight line on the axis of Zq. 

Now, let us decompose F„ into its three parts: 

pA' pB' pB 

Fn= e"'^(")d2+ / e"^(^)d;z+ / e'^^^'Mz. (138) 

J A J A' Jb' 



^"Notc that in the direction 9 — Tr/2, which is perpendicular to the axis, a.rg[h" {zo){z — zo)^] = tt — tt = 0, 
which means that h"(zo)(z — zq)^ is real and positive (i.e., it behaves like a positive parabola). Therefore, 
in this direction, zq is a local minimum. 
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B 



B' 



axis 



A! 



Figure 4: A path V from ^4 to i?, passing via zq along the axis. 



As for the first and the third terms, 



A' pB 



J +J ]dz\e''''^'^\=IJ +j j dze""(^) (139) 



< 



whose contribution is neghgible compared to e""(^"\ just hke the tails in the Laplace method. 
As for the middle integral, 

-B' rB' 



/ ^nh{z)^^^^nh{zo) / exp{n/i"(^o)(^-^o)V2}d;: 

J A' J A' 



(140) 



By changing from the complex integration variable z to the real variable a;, running from —5 
to +5, with z — zq-\- xe^^ (motion along the axis), we get exactly the Gaussian integral of 
the Laplace method, leading to 

r-B' 



expM'(..)(.-..)72}d. = e^y-^ 
where the factor e^^ is due to the change of variable {dz — e^^dx). Thus, 



(141) 



27r 



n\h"{z^)\ 



(142) 



and slightly more generally. 
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The idea of integration along the axis is that along this direction, the 'phase' of e"'**^^) is locally 
constant, and only the modulus varies. Had the integration been along another direction 
with an imaginary component j(f){z), the function e"^*^^) would have undergone 'modulation', 
i.e., it would have oscillated with a complex exponential e^^'^^^^ of a very high 'frequency' 
(proportional to n) and then e""^^°) would not have guaranteed to dictate the modulus and 
to dominate the integral. 

Now, an important comment is in order: What happens if there is more than one saddle 
point? Suppose we have two saddle points, zi and Z2- On a first thought, one may be 
concerned by the following consideration: We can construct two paths from Ato B, path Vi 
crossing Zi, and path 7^2 crossing Z2- Now, if Zi is the highest point along Vi for both i — 1 
and i = 2, then F„ is exponentially both e"'**^^!^ and e"''**^^^^ at the same time. If h{zi) ^ h{z2), 
this is a contradiction. But the following consideration shows that this cannot happen as 
long as h{z) is analytic within the region C surround by Vi U V2- Suppose conversely, that 
the scenario described above happens. Then either zi or Z2 maximize |e'*'^*^^)| along the closed 
path 7^1 U 7^2- Let us say that it is zi. We claim that then zi cannot be a saddle point, for 
the following reason: No point in the interior of C can be higher than zi, because if there 
was such a point, say, z^, then we had 

max|e"''(^)| > |e"''(^3)| > unM^i)! = ^^x \e'"''-'^\ (143) 

which contradicts the maximum modulus principle. This then means, among other things, 
that in every neighborhood of zi, all points in C are lower than Zi, including points found in a 
direction perpendicular to the direction of the axis through Zi. But this contradicts the fact 
that Zi is a saddle point: Had it been a saddle point, it would be a local maximum along the 
axis and a local minimum along the perpendicular direction. Since zi was assumed a saddle 
point, then it cannot be the highest point on Vi, which means that it doesn't dominate the 
integral. 

One might now be concerned by the thought that the integral along Vi is then dominated 
by an even higher contribution, which still seems to contradict the lower exponential order 
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of e"'*^^^) attained by the path V2- However, this is not the case. The highest point on 
the path is guaranteed to dominate the integral only if it is a saddlepoint. Consider, for 



the modulus (or attitude) is e"" everywhere. If the attitude alone had been whatever counts 
(regardless of whether it is a saddle point or not), the exponential order of (the modulus of) 
this integral would be e"". However, the true value of this integral is zero! The reason for 
this disagreement is that there is no saddle point along this path. 

What about a path V that crosses both zi and ^2? This cannot be a good path for the 
saddle point method, for the following reason: Consider two slightly perturbed versions of 
V: path V\i which is very close to T', it crosses Zi, but it makes a tiny detour that bypasses 
^2, and similarly path V2, passing via Z2, but with a small deformation near Zi. Path V2 
includes Z2 as saddle point, but it is not the highest point on the path, since V2 passes near 
2^1, which is higher. Path Vi includes Zi as saddle point, but it cannot be the highest point on 
the path because we are back to the same situation we were two paragraphs ago. Since both 
Vi and V2 are bad choices, and since they are both arbitrarily close to V, then V cannot be 
good either. 

To summarize: if we have multiple saddle points, we should find the one with the lowest 
attitude and then we have a chance to find a path through this saddlepoint (and only this 
one) along which this saddle point is dominant. 

Let us look now at a few examples. 

Example 1 - relation between D,{E) and Z{P) revisited. Assuming, without essential loss 
of generality, that the ground-state energy of the system is zero, we have seen before the 
relation Z{I3) = dEQ{E)e~^^ , which actually means that Z((3) is the Laplace transform 
of fl{E). Consequently, this means that fl{E) is the inverse Laplace transform of Z{f3), i.e.. 



where the integration in the complex plane is along the vertical line Re(/3) = 7, which is 



example, the integral F„ = 



fa+j2-n- 
a+jO 



e"'^dz. Along the vertical line from a + jO to a + j27r, 




(144) 
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chosen to the right of all singularity points of Z{(3). In the large n limit, this becomes 

fl{E) e^^^'+'^^^Mp, (145) 

which can now be assessed using the saddle point method. The derivative of the bracketed 
term at the exponent vanishes at the value of j3 that solves the equation 0'(/3) = — e, which is 
^*(e) e IR, thus we will choose 7 = /5*(e) (assuming that this is a possible choice) and thereby 
let the integration path pass through this saddle point. At ^ — /5*(e), | cxp{n[/3e + 0(/3)]}| 
has its maximum along the vertical direction, /3 = /3*(e) + jcj, —00 < u < +00 (and 
hence it dominates the integral), but since it is a saddle point, it minimizes | exp{n[/3e + 
0(/3)]}| = exp{n[/3e + 0(/3)]}, in the horizontal direction (the real hne). Thus, Q,{E) — 
exp{nminy3g]R[^e + 0(/5)]} = e^^^^\ as we have seen before. 

Example 2 - size of a type class. Here is a question which we know how to answer using 
the method of types. Among all binary sequences of length A?", how many have n I's and 



{N - n) O's? 



a;e{o,i}^ I j=i 

1 1 ( N 

a;i=0 a;jv=0 i=l 



Xi 



2^ , ^ 
27r 

i=l 



.Xi=0 
/■27r J 

= / _e^'"'^(l + e--^''^)^ 
Jo 27r 

= /'^^exp{iVb-a;a + ln(l + e-^-'^)]} 
/■^''^ d^ 

= / — exp{A^[^Q; + ln(l + e-^)]} ju — >z (146) 
Jo 27rj 
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This is an integral with a starting point A at the origin and an ending point B at 27rj. 
Here, h{z) — za + ln(l + e~^), and the saddle point, where h'{z) — 0, is on the real axis: 
Zo = In^^, where h{zo) gives the binary entropy of a, as expected. Thus, the integration 
path must be deformed to pass through this point on the real axis, and then to approach 
back the imaginary axis, so as to arrive at B. There is one serious caveat here, however: The 
points A and B are both higher than zq: While u{zo) — — Q;ln(l — a) — {1 — a) ln(l — a), at 
the edges we have u{A) = u{B) — In 2. So this is not a good saddle-point integral to work 
with. 

Two small modifications can, however, fix the problem: The first is to define the inte- 
gration interval of uu to be [— tt, tt] rather than [0, 27r] (which is, of course, legitimate), and 
then z would run from — jtt to +jTr. The second is the following: Consider again the first 
line of the expression of M„ above, but before we do anything else, let us multiply the whole 
expression (outside the summation) by e^" {9 an aribtrary real), whereas the summand will 
be multiplied by e~^^<^% which exactly cancels the factor of e^" for every non-zero term of 
this sum. We can now repeat exactly the same calculation as above (exercise), but this time 
we get: 



namely, we moved the integration path to a parallel vertical line and shifted it by the 
amount of tt to the south. Now, we have the freedom to choose 9. The obvious choice is 



path (exercise: please verify). Moreover, the vertical direction of the integration is also the 
direction of the axis of Zq (exercise: verify this too), so now everything is fine. Also, the 
second order factor of 0{l/y/n) of the saddle point integration agrees with the same factor 
that we can see from the Stirling approximation in the more refined formula. 

A slightly different look at this example is as follows. Consider the Schottky example 
and the partition function 




(147) 



to set ^ = Ini^, 



SO that we cross the saddle point Zq- Now is the highest point on the 




(148) 
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which, on the one hand, is given by X]^=o ^^'^ other hand, is given also by 

(1 + 6-^3^0)^. Thus, defining s = 6'^'°, we have Z{s) = J2n=o ^nS"", and so, Z{s) = + 
is the 2;-transform of the finite sequence {Mn}n=o- Consequently, Mn is given by the inverse 



^-transform of Zfs) 



[l + s 



I.e., 



\N-n-l 



ds 



= — /exp{A^[ln(l + s) - Q;lns]}ds 

27rj J 



(149) 



This time, the integration path is any closed path that surrounds the origin, the saddle point 
is So = ce/(l ~ ct)) so we take the path to be a circle whose radius is r = ^z^- The rest of 
the calculation is essentially the same as before, and of course, so is the result. Note that 
this is actually the very same integral as before up to a change of the integration variable 
from z to s, according to s = e~^, which maps the vertical straight line between 6 — ttj and 
9 + TTj onto a circle of radius e~^, centered at the origin. □ 

Example 3 - surface area of a sphere. Let us compute the surface area of an n-dimensional 
sphere with radius nR: 

Sn ^ I dx8[nR-y^. 

I dxe""''^*^^ ■ 6 nR — \^ ] (a > to be chosen later.) 



^naR 



, +00 1/1 

2 / ^idinR-JZixj) 



27r 



— e 



naR 



jnaR 



oo 
oo 



-oo 

+ 00 



jenR 



TT 



oo 

n/2 r+oo 



d9 
— t 

'^^ J8nR 

27r \a+j9 



TT 



n/2 
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d9 exp < n 



TT 



— oo 
n/2 pa+joo 



{a + j9)R--\n{a + j9) 



27T 



dz exp < n 



a—joo 



zR In z 

2 



(150) 



67 



So here h{z) — zR — | In 2; and the integration is along an arbitrary vertical straight line 
parametrized by a. We will choose this straight line to pass thru the saddle point Zq — 

(exercise: show that this is indeed the highest point on the path). Now, h{zQ) = ^ ln(27rei?), 
just like the differential entropy of a Gaussian RV (is this a coincidence?). □ 

Comment: In these examples, wc used an additional trick: whenever we had to deal with an 
'ugly' function like the 6 function, we presented it as an inverse transform of a 'nice' function, 
and then changed the order of integrations/summations. This idea will be repeated in the 
sequel. It is used very frequently by physicists. 

3.4 The Replica Method 

The replica method is one of the most useful tools, which originally comes from statistical 
physics, but it finds its use in a variety of other fields, with Communications and Informa- 
tion Theory included (e.g., multiuser detection). As we shall see, there are many models 
in statistical physics, where the partition function Z depends, among other things, on a 
bunch of random parameters (to model disorder), and then Z, or InZ, becomes, of course, 
a random variable as well. Further, it turns out that more often than not, the RV - In Z 
exhibits a concentration property, or in the jargon of physicists, a self-averaging property: 
in the thermodynamic limit of n — )> 00, it falls in the vicinity of its expectation ^(In Z), with 
very high probability. Therefore, the computation of the per~particle free energy (and hence 
also many other physical quantities), for a typical realization of these random parameters, 
is associated with the computation of (In Z) . The problem is that in most of the interest- 
ing cases, the exact closed form calculation of this expectation is extremely difficult if not 
altogether impossible. This is the point where the replica method enters into the picture. 

Before diving into the description of the replica method, it is important to make a certain 
digression: This is a non-rigorous, heuristic method, and it is not quite clear (yet) what are 
exactly the conditions under which it gives the correct result. Physicists tend to beheve in 
it very strongly, because in many situations it gives results that make sense, live in harmony 
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with intuition, or make good fit to experimental results and/or simulation results. The 
problem is that when there are no other means to test its validity, there is no certainty that 
it is credible and reliable. In such cases, 1 believe that the correct approach would be to refer 
to the results it provides, as a certain educated guess or as a conjecture, rather than a solid 
scientific truth. As we shall see shortly, the problematics of the rephca method is not just 
that it depends on a certain interchangeability between a limit and an integral, but more 
severely, that the procedure that it proposes, is actually not even well-defined. In spite of 
all this, since this method is so widely used, it would be inappropriate to completely ignore 
it in a course of this kind, and therefore, we will devote to the replica method at least a 
short period of time, presenting it in the general level, up to a certain point. However, we 
will not use the replica method elsewhere in this course. 

Consider then the calculation of E In Z. The problem is that Z is a sum, and it is not 
easy to say something intelligent on the logarithm of a sum of many terms, let alone the 
expectation of this log-sum. If, instead, we had to deal with integer moments of Z, EZ"^, 
life would have been much easier, because integer moments of sums, are sums of products. 
Is there a way then that we can relate moments EZ'^ to E In Zl The answer is, in principle, 
affirmative if real, rather than just integer, moments are allowed. These could be related 
via the simple relation 

, EZ"^ - 1 , In EZ"" , , 

i;inZ=lim = lim 151 

provided that the expectation operator and the limit over m can be interchanged. But we 
know how to deal only with integer moments of m. The first courageous idea of the replica 
method, at this point, is to offer the following recipe: Compute EZ'^, for positive integer 
m, and obtain an expression which is a function of m. Once this has been done, now forget 
that m is an integer, and think of it as a real variable. Finally, use the above identity, taking 
the limit of m ^ 0. 

Beyond the technicality of interchanging the expectation operator with the limit, which 
is, after all, OK in most conceivable cases, there is a more serious concern here, and this is 
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that the above procedure is not well-defined, as mentioned earlier: We derive an expression 
f{m) = EZ"^, which is originally meant for m integer only, and then 'interpolate' in between 
integers by using the same expression, in other words, we take the analytic continuation. 
Actually, the right-most side of the above identity is /'(O) where /' is the derivative of /. 
But there are infinitely many functions of a continuous variable m that pass through given 
points at integer values of m: If f{m) is such, then /(m) = fim) + g{m) is good as well, 
for every g that vanishes on the integers, for example, take g{'m) = Asm{'iirn). Nonetheless, 
/'(O) might be different from /'(O), and this is indeed the case with the example where g 
is sinusoidal. So in this step of the procedure there is some weakness, but this is simply 
ignored... 

After this introduction, let us now present the replica method on a concrete example, 
which is essentially taken from the book by Mezard and Montanari. In this example, Z = 
EL where {Ei}^l ^ are i.i.d. RV's. In the sequel, we will work with this model quite 

a lot, after we see why, when and where it is relevant. It is called the random energy model 
(REM). But for now, this is just a technical example on which we demonstrate the rephca 
method. As the replica method suggests, let's first look at the integer moments. First, what 
we have is: _ 

2" 2" m 

= J2...J2exp{-P^E,J. (152) 

The right-most side can be thought of as the partition function pertaining to a new system, 
consisting of m independent replicas (hence the name of the method) of the original system. 
Each configuration of the new system is indexed by an m-tuple i — (ii, . . . , i^), where each 



2" 




i=l 
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ia runs from 1 to 2", and the energy is Let us now rewrite Z'^ slightly differently: 

2" 2" ( m 



21=1 im=^ \ a=l 

m 2" 



exp < —(3 = j)Ej > !(•) = indicator function 

a=l j=l 



2" m 



= ^exp<^-/3 5^5^X(i„=j)^,- 

j I, j=l a=l 

2" r m 

= En^^pi-/^E^(^'^=^')^^- 

Let us now further suppose that each Ej is A/'(0, nJ^/2), as is customary in the REM, for 

reasons that we shall see later on. Then, taking expecations w.r.t. this distribution, we get: 

2" r jTt 
EZ^ = Y,El[expi~(3j2li^a^m 

= exp < — ^ — X(ia = j)I{ib = j)f using independence and Gaussianity 

i J=l I a,6=l 

= E ^^p ] E E = ^^^^^f> = 

j I a,b=l j=l 

{^2^j2 "I 
E = f ■ 
o,6=l ) 

We now define an m x m binary matrix Q, called the overlap matrix, whose entries are 
Qab = ^{ia = ib)- Note that the summand in the last expression depends on i only via 
Q. Let Nn{Q) denote the number of configurations {i} whose overlap matrix is Q. We 
have to exhaust all possible overlap matrices, which are all binary symmetric matrices with 
I's on the main diagonal. Observe that the number of such matrices is 2"*^"*"^)/^ whereas 
the number of configurations is 2""*. Thus we are dividing the exponentially large number 
of configurations into a relatively small number (independent of n) of equivalence classes, 
something that rings the bell of the method of types. Let us suppose, for now, that there is 
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some function s{Q) such that Nn{Q) — e"*^'^^ and so 

= ^ e"^(^) (153) 
Q 

with: 

9{Q)^^^Qa, + s{Q). (154) 

o,6=l 

Prom this point onward, the strategy is to use the saddle point method. Note that the 
function g{Q) is symmetric under rephca permutations: let tt be a permutation operator of 
m objects and let be the overlap matrix with entries Q^f, = Qn{a)TT{b)- Then, g{Q'^) = g{Q)- 
This property is called replica symmetry (RS), and this property is inherent to the replica 
method. In light of this, the first natural idea that comes to our mind is to postulate that 
the saddle point is symmetric too, in other words, to assume that the saddle-point Q has 
I's on its main diagonal and all other entries are taken to be the same (binary) value, call it 
qq. Now, there are only two possibilities: 

• go = and then Ar„(Q) = 2"(2" - 1) • ■ • (2" -m + 1), which implies that s{Q) = mln2, 
and then g{Q) = goiQ) = m(/3VV4 + In 2), thus {In EZ'^)/m = /3VV4 + ln2, and so 
is the limit as m — > 0. Later on, we will compare this with the result obtained from a 
more rigorous derivation. 

• qo — 1, which means that all components of i are the same, and then Nn{Q) — 2"^, 
which means that s{Q) —\n2 and so, g{Q) — gi{Q) = m^^^J^/4 + In 2. 

Now, one should check which one of these saddle points is the dominant one, depending on 
f3 and m. For m > 1, the behavior is dominated by max{go{Q) , gi{Q)} , which is gi{Q) for 
f3 > Pcim) = jA/ln 2/m, and go{Q) otherwise. For m < 1 (which is, in fact, the relevant 
case for m 0), one should look at mm{go{Q) , gi{Q)} (!), which is go{Q) in the high- 
temperature range. As it turns out, in certain regions in the /3-m plane, we must back off 
from the 'belief that dominant configurations are purely symmetric, and resort to the quest 
for dominant configurations with a lower level of symmetry. The first step, after having 
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exploited the purely symmetric case above, is called one-step replica symmetry breaking 
(IRSB), and this means some partition of the set {1,2, ... ,m} into two complementary 
subsets (say, of equal size) and postulating a saddle point Q of the following structure: 



In further steps of symmetry breaking, one may split {1,2,..., m} to a larger number of 
subsets or even introduce certain hierarchical structures. The replica method includes a 
variety of heuristic guidelines in this context. We will not delve into them any further in the 
framework of this course, but the interested student/reader can easily find more details in 
the literature, specifically, in the book by Mezard and Montanari. 




1 a = fe 

go 0, and b are in the same subset 
Qi a and b are in different subsets 
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4 Interacting Particles and Phase Transitions 

4.1 Introduction — Origins of Interactions 

As I said already in the introductory part on the analysis tools and asymptotic methods, 
until now, we have dealt almost exclusively with systems that have additive Hamiltonians, 
S{x) = '^i£{xi), which means that the particles are i.i.d. and there is no interaction: each 
particle behaves as if it was alone in the world. In Nature, of course, this is seldom really 
the case. Sometimes this is still a reasonably good approximation, but in many others the 
interactions are appreciably strong and cannot be neglected. Among the different particles 
there could be many sorts of mutual forces, e.g., mechanical, electrical, magnetic, etc. There 
could also be interactions that stem from quantum-mechanical effects: Pauli's exclusion 
principle asserts that for a certain type of particles, called Fermions (e.g., electrons), no 
quantum state can be populated by more than one particle. This gives rise to a certain 
mutal influence between particles. Another type of interaction stems from the fact that the 
particles are indistinguishable, so permutations between them are not considered as distinct 
states. We have already seen this as an example at the beginning of the previous set of 
lecture notes: In a quantum gas, as we eliminated the combinatorial factor (that counted 
indistinguishable states as distinguishable ones), we created statistical dependence, which 
physically means interactions]^ 

4.2 A Few Models That Will be Discussed in This Subsection 
Only 

The simplest forms of deviation from the purely additive Hamiltonian structure are those 

that consists, in addition to the individual energy terms {S{xi)}, also terms that depend on 

pairs, and/or triples, and/or even larger cliques of particles. In the case of purely pairwise 

"'^"'^Indeed, in the case of the boson gas, there is a weh-known effect referred to as Bose-Einstein conden- 
sation, which is actuahy a phase transition, but phase transitions can occur only in systems of interacting 
particles, as will be discussed in this set of lectures. 
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interactions, this means a structure like the following: 

n 

S{x) = '^£{xi) + ^e{xi,Xj) (156) 

where the summation over pairs can be defined over all pairs i j, or over some of the pairs, 
according to a given rule, e.g., depending on the distance between particle i and particle 
j, and according to the geometry of the system, or according to a certain graph whose 
edges connect the relevant pairs of variables (that in turn, are designated as nodes). For 
example, in a one-dimensional array (a lattice) of particles, a customary model accounts for 
interactions between neighboring pairs only, neglecting more remote ones, thus the second 
term above would be £(xj, Xj+i). A well known special case of this is that of a solid, i.e., 
a crystal lattice, where in the one-dimensional version of the model, atoms are thought of 
as a chain of masses connected by springs (see left part of Fig. [5]), i.e., an array of coupled 
harmonic oscillators. In this case, e{xi, Xj+i) = ^K{ui+i — Ui)'^, where is a constant and Ui 
is the displacement of the i-th atom from its equilibrium location, i.e., the potential energies 
of the springs. This model has an easy analytical solution (by applying a Fourier transform 
on the sequence {««}), where by "solution", we mean a closed-form, computable formula 
for the log-partition function, at least in the thermodynamic limit. In higher dimensional 




Figure 5: Elastic interaction forces between adjacent atoms in a one-dimensional lattice (left part 
of the figure) and in a two-dimensional lattice (right part). 
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arrays (or lattices), similar interactions apply, there are just more neighbors to each site, 
from the various directions (see right part of Fig. In a system where the particles are 
mobile and hence their locations vary and have no geometrical structure, like in a gas, the 
interaction terms are also potential energies pertaining to the mutual forces (see Fig. |6]), and 
these normally depend solely on the distances \\fi — fj\\. For example, in a non-ideal gas. 




-0: 




Figure 6: Mobile particles and mutual forces between them. 



i=l 



(157) 



A very simple special case is that of hard spheres (Billiard balls), without any forces, where 



oo llr,- — fA\ < 2R 







If,- - fA\ > 2R 



158) 



which expresses the simple fact that balls cannot physcially overlap. This model can (and 
indeed is) being used to obtain bounds on sphere-packing problems, which are very relevant 
to channel coding theory. This model is also solvable, but this is beyond the scope of this 
course. 

4.3 Models of Magnetic Materials — General 



Yet another example of a model, or more precisely, a very large class of models with interac- 
tions, are those of magnetic materials. These models will closely accompany our dicussions 
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from this point onward, because some of them lend themselves to mathematical formalisms 
that are analogous to those of coding problems, as we shall see. Few of these models are 
solvable, but most of them are not. For the purpose of our discussion, a magnetic material is 
one for which the important property of each particle is its magnetic moment. The magnetic 
moment is a vector proportional to the angular momentum of a revolving charged particle 
(like a rotating electron, or a current loop), or the spin, and it designates the intensity of its 
response to the net magnetic field that this particle 'feels'. This magnetic field may be the 
superposition of an externally applied magnetic field and the magnetic fields generated by 
the neighboring spins. 

Quantum mechanical considerations dictate that each spin, which will be denoted by Sj, 
is quantized - it may take only one out of finitely many values. In the simplest case to 
be adopted in our study - only two values. These will be designated by Sj = +1 ("spin 
up") and Si = —1 ("spin down"), corresponding to the same intensity, but in two opposite 
directions, one parallel to the magnetic field, and the other - antiparallel (see Fig. [7]). The 



Figure 7: Ilustration of a spin array on a square lattice. 

Hamiltonian associated with an array of spins s = (si, . . . , s„) is customarily modeled (up 
to certain constants that, among other things, accommodate for the physical units) with a 
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structure like this: 

n 

S{s) = -B - ^Si-^JijSiSj, (159) 

where B is the externally applied magnetic field and {Jij} are the coupling constants that 
designate the levels of interaction between spin pairs, and they depend on properties of 
the magnetic material and on the geometry of the system. The first term accounts for the 
contributions of potential energies of all spins due to the magnetic field, which in general, 

— * 

are given by the inner product B ■ Si, but since each Si is either parallel or antiparallel to 

— * 

B, as said, these boil down to simple products, where only the sign of each Si counts. Since 
P{s) is proportional to e~^^^^\ the spins 'prefer' to be parallel, rather than antiparallel to 
the magnetic field. The second term in the above Hamiltonian accounts for the interaction 
energy. If J^j are all positive, they also prefer to be parallel to one another (the probability 
for this is larger) , which is the case where the material is called ferromagnetic (like iron and 
nickel). If they are all negative, the material is antiferromagnetic. In the mixed case, it is 
called a spin glass. In the latter, the behavior is rather complicated, as we shall see later on. 

Of course, the above model for the Hamiltonian can (and, in fact, is being) generalized 
to include interactions formed also, by triples, quadruples, or any fixed size p (that does 
not grow with n) of spin-cliques. At this point, it is instructive to see the relation between 
spin-array models (especially, those that involve large cliques of spins) to channel codes, in 
particular, linear codes. Consider a linear code defined by a set of m partiy-check equations 
(in GF{2)), each involving the modulo-2 sum of some subset of the components of the 
codeword x. I.e., the i-th equation is: © © • • • ® xf^^ — 0, i — 1, . . . ,m. Transforming 
from Xi e {0, 1} to Si G {—1, +1} via Si — 1 — 2xi, this is equivalent to s^ts^i ■ ■ • s^i = 1. 

12 

The MAP decoder would estimate s based on the posterior 

P{s\y) = Z{y) = J] P{s)P{y\s) = P{y), (160) 

where P{s) is normally assumed uniform over the codewords (we will elaborate on this pos- 
terior later). Assuming, e.g., a BSC or a Gaussian channel P{y\s), the relevant distance 
between the codeword s = (si, . . . , Sn) and the channel output y = {yi, . . . , Hn) is propor- 
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tional to — ?/||^ = const. — 2 J2i ^iVi- Thus, P{s\y) can be thought of as a B-G distribution 
with Hamiltonian 

n m 

S{s\y) = -J^s^Vi + ^(f){s,is,i^ ■ ■ ■ SiiJ (161) 
1=1 e=i 

where J is some constant (depending on the channel parameters), the function (f){u) vanishes 
for u = 1 and becomes infinite for u ^ 1, and the partition function given by the denominator 
of P{s\y). The first term plays the analogous role to that of the contribution of the magnetic 
field in a spin system model, where each 'spin' Si 'feels' a different magnetic field proportional 
to yi, and the second term accounts for the interactions among cliques of spins. In the case 
of LDPC codes, where each parity check equation involves only a small number of bits {sj}, 
these interaction terms amount to cliques of relatively small sizeso For a general code, the 
second term is replaced by 0c (•s), which is zero for s G C and infinite otherwise. 

Another aspect of this model of a coded communication system pertains to calculations 
of mutual information and capacity. The mutual information between S and Y is, of course, 
given by 

I{S;Y) = H{Y) - H(Y\S). (162) 

The second term is easy to calculate for every additive channel - it is simply the entropy of 
the additive noise. The first term is harder to calculate: 

H{Y) = -E{\nP{Y)} = -E{\nZ{Y)}. (163) 

Thus, we are facing a problem of calculating the free energy of a spin system with random 
magnetic fields designated by the components of Y. This is the kind of calculations we 
mentioned earlier in the context of the replica method. Indeed, the replica method is used 
extensively in this context. 



^^Error correction codes can be represented by bipartite graphs with two types of nodes: variable nodes 
corresponding to the various Si and function nodes corresponding to chques. There is an edge between 
variable node i and function node j if Si is a member in clique j. Of course each Si may belong to more than 
one clique. When all cliques are of size 2, there is no need for the function nodes, as edges between nodes i 
and j simply correspond to partity check equations involving Si and Sj . 
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As we will see in the sequel, it is also customary to introduce an inverse temperature 
parameter ^, by defining 



The motivations of this will be discussed extensively later on. 

We will get back to this important class of models, as well as its many extensions, shortly. 
But before that, we discuss a very important effect that exists in some systems with strong 
interactions (both in magnetic materials and in other models) : the effect of phase transitions. 

4.4 Phase Transitions — A Qualitative Discussion 

Loosely speaking, a phase transition means an abrupt change in the collective behavior of a 
physical system, as we change gradually one of the externally controlled parameters, like the 
temperature, pressure, or magnetic field, and so on. The most common example of a phase 
transition in our everyday life is the water that we boil in the kettle when we make coffee, or 
when it turns into ice as we put it in the freezer. What exactly are these phase transitions? 
Before we refer to this question, it should be noted that there are also "phase transitions" 
in the behavior of communication systems: As the SNR passes a certain limit (for which 
capacity crosses the coding rate) , there is a sharp transition between reliable and unreliable 
communication, where the error probability (almost) 'jumps' from to 1 or vice versa. We 
also know about certain threshold effects in highly non-linear communication systems. Are 
there any relationships between these phase transitions and those of physics? We will see 
shortly that the answer is generally affirmative. 

In physics, phase transitions can occur only if the system has interactions. Consider, the 
above example of an array of spins with B = 0, and let us suppose that all Jjj > are equal. 




(164) 



where ^ controls the sharpness of the posterior distribution and 




(165) 



s 
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and thus will be denoted commonly by J. Then, 



P{s) 



(166) 



and, as mentioned earlier, this is a ferromagnetic model, where all spins 'like' to be in the 
same direction, especially when /3 and/or J is large. In other words, the interactions, in this 
case, tend to introduce order into the system. On the other hand, the second law talks about 
maximum entropy, which tends to increase the disorder. So there are two conflicting effects 
here. Which one of them prevails? 

The answer turns out to depend on temperature. Recall that in the canonical ensemble, 
equilibrium is attained at the point of minimum free energy / = e — Ts(e). Now, T plays the 
role of a weighting factor for the entropy. At low temperatures, the weight of the second term 
of / is small, and minimiizing / is approximately (and for T = 0, this is exact) equivalent to 
minimizing e, which is obtained by states with a high level of order, as £{s) = —J'Yli{ij) ^i^j, 
in this example. As T grows, however, the weight of the term —Ts{e) increases, and min/, 
becomes more and more equivalent to maxs(e), which is achieved by states with a high 
level of disorder (see Fig. |8]). Thus, the order-disorder characteristics depend primarily 



Figure 8: Qualitative graphs of /(e) at various temperatures. The minimizing e increases with T. 

on temperature. It turns out that for some magnetic systems of this kind, this transition 
between order and disorder may be abrupt, in which case, we call it a phase transition. 
At a certain critical temperature, called the Curie temperature, there is a sudden transition 
between order and disorder. In the ordered phase, a considerable fraction of the spins align in 




/ 



T = 
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the same direction, which means that the system is spontaneously magnetized (even without 
an external magnetic field), whereas in the disordered phase, about half of the spins are in 

either direction, and then the net magnetization vanishes. This happens if the interactions, 
or more precisely, their dimension in some sense, is strong enough. 

What is the mathematical significance of a phase transition? If we look at the partition 
function, Z{I3), which is the key to all physical quantities of interest, then for every finite 
n, this is simply the sum of a bunch of exponentials in (3 and therefore it is continuous 
and differentiable as many times as we want. So what kind of abrupt changes could there 
possibly be in the behavior of this function? 

It turns out that while this is true for all finite n, it is no longer necesarily true if we 
look at the thcrmodynamical limit, i.e., if we look at the behavior of = lim^^oo ^J^^^ ■ 
While must be continuous for all /3 > (since it is convex), it need not necessarily have 
continuous derivatives. Thus, a phase transition, if exists, is fundamentally an asymptotic 
property, it may exist in the thcrmodynamical limit only. While a physical system is, after 
all finite, it is nevertheless well approximated by the thcrmodynamical limit when it is very 
large. By the same token, if we look at the analogy with a coded communication system: for 
any finite block-length n, the error probability is a 'nice' and smooth function of the SNR, 
but in the limit of large n, it behaves like a step function that jumps between and 1 at the 
critical SNR. We will see that the two things are related. 

Back to the physical aspects, the above discussion explains also why a system without 
interactions, where all {xi} are i.i.d., cannot have phase transitions. In this case, Zn{f3) — 
[Zi(^)]", and so, (f){/3) — In Zi{/3), which is always a 'nice' function without any irregularities. 
For a phase transition to occur, the particles must behave in some collective manner, which 
is the case only if interactions take place. 

There is a distinction between two types of phase transitions: 

• If 4>{f3) has a discontinuous first order derivative, then this is called a first order phase 
transition. 



82 



• If has a continuous first order derivative, but a discontinuous second order deriva- 
tive then this is called a second order phase transition, or a continuous phase transition. 

We can talk, of course, about phase transitions w.r.t. additional parameters other than 
temperature. In the above magnetic example, if we introduce back the magnetic field B 
into the picture, then Z, and hence also 0, become functions of B too. If we then look at 
derivative of 



0(/3,i?)= hmi^^^M)^ linilln 



exp < ^ Si + /3 J ^ SiSj 

S I i=l {i,3) 



(167) 



w.r.t. the product (/S-B), which multiplies the magnetization, X^jSi, at the exponent, this 
would give exactly the average magnetization per spin 



m{P,B) = (^Y.s}j, 



:i68i 



and this quantity might not always be continuous. Indeed, as I mentioned earlier, below 
the Curie temperature there might be a spontaneous magnetization. If S J, 0, then this 
magnetization is positive, and if i? 0, it is negative, so there is a discontinuity at i? = 0. 
We will see this more concretely later on. We next discuss a few solvable models of spin 
arrays, with and without phase transitions. 

4.5 The One— Dimensional Ising Model 

According to this model, 

n n 

£:(s) = -Sj]si- (169) 

i=l i=l 

with the periodic boundary condition Sn+i — Si. Thus, 

{n n \ 

j3B Si + /3J^^ SjSj+i > Note: the kind of sums encountered in Markov chains 
i=\ 1=1 ) 

f 71 n ^ 

= ^exp J /i^Si + is:^SiSi+i I h = PB, K = pj 



1=1 i=l 

n 



— exp I — ^^(sj + Si+i) + K SiSi+i ^ (just to symmetrize the expression) 



, 2 

S I i=l i=l 
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Consider now the 2x2 matrix P whose entries are exp{|(s + s') + Kss'}, s,s e {— 1,+1}, 
i.e., 

/ pK+h .-K \ 

Also, Si — +1 will be represented by the column vector cTj = (1,0)-^ and Sj = —1 will be 
represented by CTj = (0, 1)^. Thus, 

= ^alP ■ I ■ P ■ I ■ ■ ■ I ■ P(7i 

= tr{P"} 

= A^ + A^ (171) 
where Ai and A2 are the eigenvalues of P, which are 



Ai_2 = cosh(/i) ± -^6-2^ + smh^{h). (172) 
Letting Ai denote the larger (the dominant) eigenvalue, i.e., 

Ai = cosh(/i) + ^e-2^^ + e2-^ sinh^(/i), (173) 

then clearly, 

0(/i, K) = lim ^ = In Ai. (174) 

n— ^-oo 77, 

The average magnetization is 
M{h,K) = 

\i=i 



ainz(/t,is:) 

dh 



(175) 
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and so, the per-spin magnetization is: 



n^-oo n oh 



or, returning to the original parametrization: 



, ^ sinh(/3i?) , , 

m(/3, 5) = , (177) 

^e-4/3-^ + sinh2(/35) 

For /3 > and i? > this is a nice function, and so, there is are no phase transitions and 
no spontaneous magnetization at any finite temperatureo However, at the absolute zero 
(/3 — )■ oo), we get 

lim lim m(/3,5) = +1; lim lim m(/3, 5) = -1, (178) 

BIO jS-^oo _BtO /3— 5>oo 

thus m is discontinuous w.r.t. i? at /3 — )■ oo, which means that there is a phase transition at 
T = 0. In other words, the Curie temperature is Tc = 0. 

We see then that one-dimensional Ising model is easy to handle, but it is not very 
interesting in the sense that there is actually no phase transition. The extension to the 
two-dimensional Ising model on the square lattice is surprisingly more difficult, but it is still 
solvable, albeit without a magnetic field. It was first solved by Onsager in 1944, who has 
shown that it exhibits a phase transition with Curie temperture given by 

Tc = , (179) 

fcln(v^+l) ^ ' 

where k is Boltzmann's constant. For lattice dimension > 3, the problem is still open. 

It turns out then that whatever counts for the existence of phase transitions, is not the 
intensity of the interactions (designated by the magnitude of J), but rather the "dimension- 
ality" of the structure of the pairwise interactions. If we denote by the number of £-th 
order neighbors of every given site, namely, the number of sites that can be reached within 
£ steps from the given site, then whatever counts is how fast does the sequence {ng} grow. 



-"^^Note, in particular, that for J = (i.i.d. spins) we get paramagnetic characteristics to(/3, B) ~ tanh(/3i?), 
in agreement with the resuh pointed out in the example of two-level systems, in one of our earlier discussions. 



85 



or more precisely, what is the value oi d — lim£_^oo i 



Inne, which is exactly the ordinary 



dimensionality for hypercubic lattices. Loosely speaking, this dimension must be sufficiently 
large for a phase transition to exist. 

To demonstrate this point, we next discuss an extreme case of a model where this dimen- 
sionality is actually infinite. In this model "everybody is a neighbor of everybody else" and 
to the same extent, so it definitely has the highest connectivity possible. This is not quite a 
physically realistic model, but the nice thing about it is that it is easy to solve and that it 
exhibits a phase transition that is fairly similar to those that exist in real systems. It is also 
intimately related to a very popular approximation method in statistical mechanics, called 
the mean field approximation. Hence it is sometimes called the mean field model. It is also 
known as the Curie-Weiss model or the infinite range model. 

Finally, I should comment that there are other "infinite-dimensional" Ising models, like 
the one defined on the Bethe lattice (an infinite tree without a root and without leaves), 
which is also easily solvable (by recursion) and it also exhibits phase transitions (see Baxter's 
book), but we will not discuss it here. 



Here, all pairs {{si,Sj)} "talk to each other" with the same "voice intensity", J/(2n), and 
without any geometry. The 1/n factor here is responsible for keeping the energy of the 
system extensive (linear in n), as the number of interaction terms is quadratic in n. The 
factor 1/2 compensates for the fact that the summation over i ^ j counts each pair twice. 
The first observation is the trivial fact that 



4.6 The Curie-Weiss Model 



According to the Curie-Weiss (C-W) model. 
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where the second equahty holds since = 1. It follows then, that our Hamiltonian is, upto 
a(n immaterial) constant, equivalent to 



J 



1=1 



-n 



,1=1 



'182) 



thus S{s) depends on s only via the magnetization m(s) = jiYli^i- This fact makes the 
C-W model very easy to handle similarly as in the method of types: 

J 



B ■ m(s) + —m^{s) 



s ^ 



m=—l 
+1 



gnh2{(l+m)/2) _ gn/3{Bm+JmV2) 



m=—l 



exp < n ■ max 

|m|<l 



, , 1 + m\ ^ ^ Bm'^J 



and so. 



B) = max 

|r?i|<l 



h2 — — + (3Bm + 



2 J 2 
The maximum is found by equating the derivative to zero, i.e.. 



1 , / 1 — m 
= - In 

2 \l + m 



(3B + (3Jm = - tanh"^ (m) + PB + (3Jm 



;i83) 



;i84) 



or equivalently, the maximizing (and hence the dominant) m is a solution m* to the equa- 
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tion 

m = tanh(/3i? + (3Jm). 
Consider first the case 5 = 0, where the equation boils down to 

m = tanh(/3 Jm). 



;i85) 



It is instructive to look at this equation graphically. Referring to Fig. [9l we have to make a 
distinction between two cases: If /3J < 1, namely, T > = J/k, the slope of the function 



^^Once again, for J = 0, we are back to non-interacting spins and then this equation gives the paramagnetic 
behavior m = tanh(/3i3). 
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y = tanh(/3Jm) at the origin, /3J, is smaller than the slope of the linear function y = m, 
which is 1, thus these two graphs intersect only at the origin. It is easy to check that in this 
case, the second derivative of ip{m) = /i2((l + m)/2) + pjrn? /2 at m = is negative, and 
therefore it is indeed the maximum (see Fig.[TOl left part). Thus, the dominant magnetization 
is m* = 0, which means disorder and hence no spontaneous magnetization for T > T^. On 



Figure 9: Graphical solutions of equation m = tanh(/3 Jm): The left part corresponds to the case 
/3J < 1, where there is one solution only, m* = 0. The right part corresponds to the case /3J > 1, 
where in addition to the zero solution, there are two non-zero solutions m* = zizniQ. 

the other hand, when /3J > 1, which means temperatures lower than Tc, the initial slope of 
the tanh function is larger than that of the linear function, but since the tanh cannot take 
values outside the interval (—1, +1), the two functions must intersect also at two additional, 
symmetric, non-zero points, which we denote by +mo and —mo (see Fig. [HI right part). In 
this case, it can readily be shown that the second derivative of iplm) is positive at the origin 
(i.e., there is a local minimum at m = 0) and negative at m = ±mo, which means that there 
are maxima at these two points (see Fig. [TUl right part). Thus, the dominant magnetizations 
are ±mo, each capturing about half of the probability. 

Consider now the case /3J > 1, where the magnetic field B is brought back into the 
picture. This will break the symmetry of the right graph of Fig. [10] and the corresponding 
graphs of iplm) would be as in Fig. [TH where now the higher local maximum (which is also 
the global one) is at mo(i?) whose sign is as that of B. But as i? — t- 0, mo(-B) niQ of Fig. 
[TOl Thus, we see the spontaneous magnetization here. Even after removing the magnetic 
field, the system remains magnetized to the level of mg, depending on the direction (the 
sign) of B before its removal. Obviously, the magnetization m(/3, B) has a discontinuity at 




y = m 



y = m 
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Figure 10: The function ^p{m) = /i2((l + + pjm? /2 lias a unique maximum at m = when 

f3J < 1 (left graph) and two local maxima at ±mo, in addition to a local minimum at m = 0, when 
(3 J > 1 (right graph). 



tp{m) 





Figure 11: The case f3J > 1 with a magnetic field B. The left graph corresponds to i? < and 
the right graph - to i? > 0. 



5 = for T < Tc, which is a first order phase transition w.r.t. B (see Fig. [T2|) . We note 
that the point T = Tc is the boundary between the region of existence and the region of 
non-existence of a phase transition w.r.t. B. Such a point is called a critical point. The 
phase transition w.r.t. /3 is of the second order. 

Finally, we should mention here an alternative technique that can be used to analyze 
this model, which is useful in many other contexts as well. It is based on the idea of using 
a transform integral, in this case, the Hubbard-Stratonovich transform, and then the saddle 
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m{P,B) 




Figure 12: Magnetization vs. magnetic field: For T < Tc there is spontaneous magnetization: 
lim^j^o iTT'iPj -B) = +mo and lims^o iT^iP^ B) = —mo, and so there is a discontinuity at 5 = 0. 



point method. Specifically, we have the following chain of equalities: 

2' 



Z{h,K 



S i=l \i=l 

= I^exp|/i|]s,|-exp|^(5^ 



h = {3B, K = /3J 



i=l 



n 



i=l 



27rK 



dz exp 



nz 
'2K 



s=-l 



n 



2tiK 



[ d^e-"^'/('^)[2cosh(/i + z)]' 



= 2" • W / d^exp{n[lncosh(/i + z) - zy{2K)]} 
V 27rA 

Using the the saddle point method (or the Laplace method), this integral is dominated by 
the maximum of the function in the square brackets at the exponent of the integrand, or 
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equivalently, the minimum of the function 

= _ -lncosh(/i + ^). (186) 

by equating its derivative to zero, we get the very same equation as m = tanh(/3S + f3Jm) 
by setting z — j3Jm. The function 7(2;) is different from the function that we maximized 
earher, but the extremum is the same. This function is called the Landau free energy. 

4.7 Spin Glass Models With Random Parameters and Random 
Code Ensembles 

So far we discussed only models where the non-zero coupling coefficients, J — {Jij} are equal, 

thus they arc cither all positive (ferromagnetic models) or all negative (antiferromagnetic 
models). As mentioned earlier, there are also models where the signs of these coefficients are 
mixed, which are called spin glass models. 

Spin glass models have a much more complicated and more interesting behavior than 
ferromagnets, because there might be metastable states due to the fact that not necessarily 
all spin pairs {(sj, Sj)} can be in their preferred mutual polarization. It might be the case that 
some of these pairs are "frustrated." In order to model situations of amorphism and disorder 
in such systems, it is customary to model the coupling coeffcients as random variables. 

Some models allow, in addition to the random coupling coefficients, also random local 
fields, i.e., the term — -B Si in the Hamiltonian, is replaced by — ^jSjSj, where {Bi} 
are random variables, similarly as in the representation of P{s\y) pertaining to a coded 
communicaion system, as discussed earlier, where {y^} play the role of local magnetic fields. 
The difference, however, is that here the {-Bj} are normally assumed i.i.d., whereas in the 
communication system model P{y) exhibits memory (even if the channel is memoryless) 
due to memory in P{s). Another difference is that in the physics model, the distribution 
of {Bi\ is assumed to be independent of temperature, whereas in coding, if we introduce 
a temperature parameter by exponentiating (i.e., Pp{s\y) oc P^{s)P^{y\s)), the induced 
marginal of y will depend on ^. 
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In the following discussion, let us refer to the case where only the coupling coefficients J 
are random variables (similar things can be said in the more general case, discussed in the 
last paragraph). This model with random parameters means that there are now two levels 
of randomness: 

• Randomness of the coupling coefficients J . 

• Randomness of the spin configuration s given J, according to the Boltzmann distri- 
bution, i.e.. 



exp |/3 




} 


Z{(5,B\J) 



P{s\J) = ^, ^. (187) 



However, these two sets of RV's have a rather different stature. The underlying setting is 
normally such that J is considered to be randomly drawn once and for all, and then remain 
fixed, whereas s keeps varying all the time (according to the dynamics of the system). At any 
rate, the time scale along which s varies is much smaller than that of J. Another difference 
is that J is normally not assumed to depend on temperature, whereas s, of course, does. 
In the terminlogy of physicists, s is considered an annealed RV, whereas J is considered a 
quenched RV. Accordingly, there is a corresponding distinction between annealed averages 
and quenched averages. 

Actually, there is (or, more precisely, should be) a parallel distinction when we consider 
ensembles of randomly chosen codes in Information Theory. When we talk about random 
coding, we normally think of the randomly chosen code as being drawn once and for all, 
we don't reselect it after each transmission (unless there are security reasons to do so), 
and so, a random code should be thought of us a quenched entity, whereas the source(s) 
and channel(s) are more naturally thought of as annealed entities. Nonetheless, this is not 
what we usually do in Information Theory. We normally take double expectations of some 
performance measure w.r.t. both source/channel and the randomness of the code, on the 
same footing]^ We will elaborate on this point later on. 



There are few exceptions to this rule, e.g., a paper by Barg and Forney, IEEE Trans, on IT, Sept. 2002, 
and several follow-ups. 
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Returning to spin glass models, let's see what is exactly the difference between the 
quenched averaging and the annealed one. If we examine, for instance, the free energy, or 
the log-partition function, In Z{(3\J), this is now a RV, of course, because it depends on the 
random J. If we denote by (■) j the expectation w.r.t. the randomness of J, then quenched 
averaging means {In Z{f5\ J)) j (with the motivation of the self-averaging property of the 
RV lnZ{/3\J) in many cases), whereas annealed averaging means ln{Z{/3\J)) j. Normally, 
the relevant average is the quenched one, but it is typically also much harder to calculate 
(and it is customary to apply the replica method then). Clearly, the annealed average is 
never smaller than the quenched one because of Jensen's inequality, but they sometimes 
coincide at high temperatures. The difference between them is that in quenched averaging, 
the dominant realizations of J are the typical ones, whereas in annealed averaging, this is 
not necessarily the case. This follows from the following sketchy consideration. As for the 
annealed average, we have: 

{z{i3\j) = Y.nj)mj) 
J 

^ ^ g-"S(o=) . g"" (assuming exponential probabihties) 

a 

_ gnmaxa[a-E(a)] (188) 

which means that the annealed average is dominated by realizations of the system with 

^^^-^^^^a* ^SiYgmeix\a-E(a)], (189) 
n a 

which may differ from the typical value of a, which is 

a = 0(/3) = lim - {In Z{p\J)) . (190) 

n— >oo n 

On the other hand, when it comes to quenched averaging, the RV In Z{/3\J) behaves hnearly 
in n, and concentrates strongly around the typical value n0(^), whereas other values are 
weighted by (exponentially) decaying probabilities. 
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In the coded communication setting, there is a strong parallehsm. Here, there is a 
distinction between the exponent of the average error probabihty. In EPe{C) (annealed) and 

the average exponent of the error probabihty ElnPe{C) (quenched), where Pe{C) is the error 
probabihty of a randomly selected code C. Very similar things can be said here too. 

The literature on spin glasses includes many models for the randomness of the coupling 
coefficients. We end this part by listing just a few. 

• The Edwards-Anderson (E-A) model, where {Jij} are non-zero for nearest-neighbor 
pairs only (e.g., j = i ± 1 in one-dimensional model). According to this model, these 
Ji/s are i.i.d. RV's, which are normally modeled to have a zero-mean Gaussian pdf, 
or binary symmetric with levels ±Jo. It is customary to work with a zero-mean 
distribution if we have a pure spin glass in mind. If the mean is nonzero, the model 
has either a ferromangetic or an anti-ferromagnetic bias, according to the sign of the 
mean. 

• The Sherrington-Kirkpatrick (S-K) model, which is similar to the E-A model, except 
that the support of {Jij} is extended to include all n{n — l)/2 pairs, and not only 
nearest-neighbor pairs. This can be thought of as a stochastic version of the C-W 
model in the sense that here too, there is no geometry, and every spin 'talks' to every 
other spin to the same extent, but here the coefficients are random, as said. 

• The p-spin model, which is similar to the S-K model, but now the interaction term 
consists, not only of pairs, but also triples, quadraples, and so on, up to cliques of size 
p, i.e., products Si^Sjj • • • Sip, where (ii, . . . , ip) exhaust all possible subsets of p spins out 
of n. Each such term has a Gaussian coefficient Ji^,...,ip with an appropriate variance. 

Considering the p-spin model, it turns out that if we look at the extreme case of p — )> oo 
(taken after the thermodynamic limit n oo), the resulting behavior turns out to be 
extremely erratic: all energy levels {£{s)} se{-i,+i}^ become i.i.d. Gaussian RV's. This is, of 
course, a toy model, which has very little to do with reality (if any), but it is surprisingly 
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interesting and easy to work with. It is called the random energy model (REM). We have 
already mentioned it as an example on which we demonstrated the replica method. We are 

next going to talk about it extensively because it turns out to be very relevant for random 
coding models. 
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5 The Random Energy Model and Random Coding 
5.1 The REM in the Absence of a Magnetic Field 



The REM was proposed by the French physicist Bernard Derrida in the early eighties of the 
previous century in a series of papers: 

1. B. Derrida, "Random-energy model: limit of a family of disordered models," Phys. 
Rev. Lett, vol. 45, no. 2, pp. 79-82, July 1980. 

2. B. Derrida, "The random energy model," Physics Reports (Review Section of Physics 
Letters), vol. 67, no. 1, pp. 29-35, 1980. 

3. B. Derrida, "Random-energy model: an exactly solvable model for disordered sys- 
tems," Phys. Rev. 5, vol. 24, no. 5, pp. 2613-2626, September 1981. 

Derrida showed in one of his papers that, since the correlations between the random energies 
of two configurations, s and s' in the p-spin model are given by 



and since I^X^ILi^^^il ^ these correlations vanish as p — >■ oo. This has motivated him 
to propose a model according to which the configurational energies {S{s)}, in the absence 
of a magnetic field, are simply i.i.d. zero-mean Gaussian RV's with a variance that grows 
linearly with n (again, for reasons of extensivity) . More concretely, this variance is taken to 
be nJ^/2, where J is a constant parameter. This means that we forget that the spin array 
has any structure of the kind that we have seen before, and we simply randomly draw an 
independent RV 6{s) ~ A/'(0,nJ^/2) (and other distributions are also possible) for every 
configuration s. Thus, the partition function Z{I3) = e^f^^(^) is a random variable as 
well, of course. 

This is a toy model that does not describe faithfully any reahstic physical system, but 
we will devote to it some considerable time, for several reasons: 




(191) 
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• It is simple and easy to analyze. 

• In spite of its simplicity, it is rich enough to exhibit phase transitions, and therefore it 
is interesting. 

• Last but not least, it will prove very relevant to the analogy with coded communication 
systems with randomly selected codes. 

As we shall see quite shortly, there is an intimate relationship between phase transitions 
of the REM and phase transitions in the behavior of coded communication systems, most 
notably, transitions between reliable and unreliable communication, but others as well. 

What is the basic idea that stands behind the analysis of the REM? As said, 

Z(/3) = ^e-'^^(*) (192) 
s 

where £{s) ~ A/'(0, nJ^/2) are i.i.d. Consider the density of states fi{E), which is now a RV: 
fl{E)dE is the number of configurations {s} whose randomly selected energy S{s) happens 
to fall between E and E + dE, and of course, 

/• + 00 
dEfl{E)e-^^. (193) 
-oo 

How does the RV fl{E)dE behave like? First, observe that, ignoring non-exponential factors: 
Pt{E < S{s) <E + dE}^ f{E)dE = c-^'/^'^'^'ME, (194) 

and so, 

{n{E)dE) = 2" • e-^'/^"''') = exp |n 

We have reached the pivotal point behind the analysis of the REM, which is based on a 
fundamental principle that goes far beyond the analysis of the first moment of Q{E)dE. In 
fact, this principle is frequently used in random coding arguments in IT: 

Suppose that we have e"^ (^4 > 0, independent of n) independent events {Si}, each one 
with probabihty Pr{£j} = e~"^ {B > 0, independent of n). What is the probabihty that 
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at least one of the 8iS would occur? Intuitively, we expect that in order to see at least one 
or a few successes, the number of experiments should be at least about 1/Pr{£^j} = e"^. If 
A > B then this is the case. On the other hand, for A < B, the number of trials is probably 
insufficient for seeing even one success. Indeed, a more rigorous argument gives: 

= i-Pr Q£:,^ 



gln(l-e-"«) 



nA 



= i_(i_e-ns)' 

= 1 

= 1 -exp{e"^ln(l -e""^)} 

^ 1 - exp{-e"^e-"^} 

= 1 -exp{-e"(^-^)} 

^ A<B (^96) 

BTW, the 2nd line could have been shown also by the union bound, as Ylii^'^i^i} = 
gjiAg-nS Exercise: What happens when A = Bl 

Now, to another question: For A > B, how many of the £j's would occur in a typical 
realization of this set of experiments? The number f2„ of 'successes' is given by X]i=i ^{^i}^ 
namely, it is the sum of e""^ i.i.d. binary RV's whose expectation is = e"^^~^\ 

Therefore, its probability distribution concentrates very rapidly around its mean. In fact, 
the events > e"(-4-B+e)j (^^ > independent of n) and < e"(-4-B-e)| g^j-e large 
deviations events whose probabilities decay exponentially in the number of experiments, 
e""^, i.e., double-exponentially (!) in Thus, for A > B, the number of successes is 

"almost deterministically" about e"^'^^^\ 

Now, back to the REM: For E whose absolute value is less than 

Eo = njVh[2 (197) 



16 



This will be shown rigorously later on. 
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the exponential increase rate, A — In 2, of the number 2" 



,n In 2 



of configurations, 



the number of independent trials in randomly drawing energies {S{s)}, is faster than the 
exponential decay rate of the probability, e""'!'^/^"'^)!^) = e""*^^/'^^^ (i.e., B = (e/J)^) that S{s) 
would happen to fall around E. In other words, the number of these trials is way larger than 
one over this probability and in view of the earlier discussion, the probability that 



n{E)dE = ^I{E < S{s) <E + dE}. 



(198) 



would deviate from its mean = exp{n[ln 2 — (E/ (nJ))^]}, by a multiplicative factor that falls 
out of the interval [e~"^, e+"^], decays double-exponentially with n. In other words, we argue 
that for —Eq < E < +Eo, the event 



exp < n 



ln2 - 



E 
nJ 



< n{E)dE < e+"' • exp <^ n 



ln2 - 



E 
nJ 



(199) 



happens with probability that tends to unity in a double-exponential rate. As discussed, 
—Eq < E < +Eq is exactly the condition for the expression in the square brackets at the 
exponent [In 2 — (^)^] to be positive, thus fl{E)dE is exponentially large. On the other 
hand, if l-B] > Eq, the number of trials 2" is way smaller than one over the probability of 
falling around E, and so, most of the chances are that we will see no configurations at all 
with energy about E. In other words, for these large values of \E\, Q{E) = for typical 
realizations of the REM. It follows then that for such a typical realization, 

-Eq 



m 



{dE ■ Q{E)) e'^^ 

■En 



Eo 
+Eo 

■Eo 
+Eo 



dE ■ exp < n 



/+Ei 
-Eo 



dE ■ exp < n 



ln2 - 



ln2 - 



E 
nJ 

E 
nJ 



n 



de ■ exp < n 



eo 



exp < n ■ max 

\e\<eo 



In 2 



-/3E 

E 
n 



In 2 



1 A -B A -Eo /—— 

> e = — , eo = — = Jvln2, 
) n n 

by Laplace integration 
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The maximization problem at the exponent is very simple: it is that of a quadratic function 
across an interval. The solution is of either one of two types, depending on whether the 
maximum is attained at a zero-derivative internal point in (— eo, +eo) or at an edgepoint. 
The choice between the two depends on {3. Specifically, we obtain the following: 



= lim 



lnZ(/3) 



n— >oo n 



ln2 + 



/3</3c 
/3jVh^ /3 > /3c 



(200) 



where Pc = jVln 2. What we see here is a phase transition. The function 0(/3) changes its 
behavior abruptly at /3 = Pc, from being quadratic in /3 to being linear in /3 (see also Fig. 
[T3l right part). The function is continuous (as always), and so is its first derivative, but 
the second derivative is not. Thus, it is a second order phase transition. Note that in the 
quadratic range, this expression is precisely the same as we got using the replica method, 
when we hypothesized that the dominant configuration is fully symmetric and is given by 



Q 



Thus, the replica symmetric solution indeed gives the correct result in the high 



temperature regime, but the low temperature regime seems to require symmetry breaking. 
Thus, the condition i? > ln2 — h2{S) is equivalent to 

m 

S(£) = ln2-(^) 





Figure 13: The entropy function and the normalized log-partition function of the REM. 



What is the significance of each one of these phases? Let's begin with the second line 



of the above expression of (j){/3), which is 0(/3) = /3 JVln 2 = /3eo for /3 > /3c- What is the 
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meaning of linear dependency of (j) in /3? Recall that the entropy S is given by 

S(/3) = 0(/3)-/3-0'(/3), 

which in the case where is linear, simply vanishes. Zero entropy means that the partition 
function is dominated by a subexponential number of ground-state configurations (with per- 
particle energy about eo), just like when it is frozen (see also Fig. [131 left part: S(— cq) = 0). 
This is why we will refer to this phase as the frozen phase or the glassy phas^^ In the high- 
temperature range, on the other hand, the entropy is strictly positive and the dominant 
per-particle energy level is e* = — |/3J^, which is the point of zero-derivative of the function 
[In 2 — (e/J)^ — (3e]. Here the partition is dominated by exponentially many (exercise: what 
is the exponent?) configurations whose energy is E* = ne* = — ^/3J^. As we shall see later 
on, in this range the behavior of the system is essentially paramagnetic (like in a system of 
i.i.d. spins), and so it is called the paramagnetic phase. 

We therefore observe that the type of phase transition here is different than in the Curie- 
Weiss model. We are not talking here about spontaneous magnetization transition, but 
rather on a glass transition. In fact, we will not see here a spontaneous magnetization even 
if we add a magnetic field (time permits, this will be seen later on). 

From 0(/3), one can go ahead and calculate other physical quantities, but we will not do 
this now. As a final note in this context, I wish to emphasize that since the calculation of Z 
was carried out for the typical realizations of the quenched RV's {S{s)}, we have actually 
calculated the quenched average of lim„(lnZ)/n. As for the annealed average, we have 



lim 

n— >oo 



Hzm 



n 



lim — In 

n— s>oo n 



{n{E)de)e 



IJTR 



lim — In 

n— >oo n 



exp < n 



TR 



l„2-l- 



max 

eelR 



i„2-ii)= 



Pi 



Laplace integration 



In 2 + 



(201) 



^^In this phase, the system behaves hke a glass: on the one hand, it is frozen (so it consohdates), but on 
the other hand, it remains disordered and amorphous, hke a hquid. 
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which is the paramagnetic expression, without any phase transition since the maximization 
over e is not constrained. 



5.2 The Random Code Ensemble and its Relation to the REM 

Let us now see how does the REM relate to random code ensembles. The discussion in 
this part is based on Mczard and Montanari's book, as well as on the paper: N. Merhav, 
"Relations between random coding exponents and the statistical physics of random codes," 
IEEE Trans. Inform. Theory, vol. 55, no. 1, pp. 83-92, January 2009. Another relevant 
paper is: A. Barg and G. D. Forney, Jr., "Random codes: minimum distances and error 
exponents," IEEE Trans. Inform. Theory, vol. 48, no. 9, pp. 2568-2573, September 2002. 

Consider a DMC, P{y\x) — Y\4=iP{yi\^i)i by an input n-vector that belongs to a 
codebook C — {xi, X2,..., Xm}, M — e"^, with uniform priors, where R is the coding rate 
in nats per channel use. The induced posterior, for x & C, is then: 

^"'"'^ " Ex'.cPiy\^') 

g-in[i/p(2/|a;)] 

" Ea;'ece"'"f'/''^^'^'^^' ^^^^^ 
Here, the second line is deliberately written in a form that resembles the Boltzmann dis- 
tribution, which naturally suggests to consider, more generally, the posterior distribution 
parametrized by ^, that is 



Y.x'ecPny\^') 

g-r^iu[i/p(?/'ic)] 



my) 

There are a few motivations for introducing the temperature parameter: 

• It allows a degree of freedom in case there is some uncertainty regarding the channel 
noise level (small /3 corresponds to high noise level). 
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• It is inspired by the ideas behind simulated anneahng techniques: by samphng from 
while gradually increasing /3 (cooling the system), the minima of the energy function 
(ground states) can be found. 

• By applying symbolwise maximum a-posteriori (MAP) decoding, i.e., decoding the 
£-th symbol of a; as argmax^ P/3(a;£ = a|t/), where 

Pp{x,^a\y)^ PMv)^ (203) 

CCeC: Xi=a 

we obtain a family oi finite-temperature decoders (originally proposed by Rujan in 1993) 
parametrized by ^, where /3 = 1 corresponds to minimum symbol error probability 

(with respect to the real underlying channel P{y\x)) and (5 ^ oo corresponds to 
minimum block error probability. 

• This is one of our main motivations: the corresponding partition function, Z{f3\y), 
namely, the sum of (conditional) probabilities raised to some power /3, is an expression 
frequently encountered in Renyi information measures as well as in the analysis of 
random coding exponents using Gallager's techniques. Since the partition function 
plays a key role in statistical mechanics, as many physical quantities can be derived 
from it, then it is natural to ask if it can also be used to gain some insights regarding 
the behavior of random codes at various temperatures and coding rates. 

For the sake of simplicity, let us suppose further now that we are dealing with the binary 
symmetric channel (BSC) with crossover probability p, and so, 

P{y\x) = p'^(^.2/)(l - py-di^'V) = (1 - p)"e-^'^(^'2/), (204) 

where J = In and d{x, y) is the Hamming distance. Thus, the partition function can be 
presented as follows: 

Z{(3\y) = (1 -pY^Ye-^^^^^^y\ (205) 
xec 

Now consider the fact that the codebook C is selected at random: Every codeword is ran- 
domly chosen independently of all other codewords. At this point, the analogy to the REM, 
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and hence also its relevance, become apparent: If each codeword is selected independently, 
then the 'energies' {Jd{x,y)} pertaining to the partition function 

Z{/3\y) = (1 ^p)PnJ2e-^-"'^'^'y\ (206) 
xec 

(or, in the case of a more general channel, the energies {— ln[l/P(^/|£c)]} pertaining to 
the partition function Z{f3\y) = Exec e"'^'"'^^^^^'^^^' 

are i.i.d. random variables for all 
codewords in C, with the exception of the codeword xq that was actually transmitted and 
generated Since we have seen phase transitions in the REM, it is conceivable to expect 
them also in the statistical physics of the random code ensemble, and indeed we will see 
them shortly. 

Further, we assume that each symbol of each codeword is drawn by fair coin tossing, i.e., 
independently and with equal probabilities for '0' and '1'. As said, we have to distinguish 
now between the contribution of the correct codeword Xq, which is 

Z,{{3\y) = (1 - j9)^"e--'^(^«'^) (207) 

and the contribution of all other (incorrect) codewords: 

Ze(/3|t/) = (l-p)^" e''''^^'y\ (208) 

xec\{Xo} 

Concerning the former, things are very simple: Typically, the channel flips about np bits out 
the n transmissions, which means that with high probability, d{xQ, y) is about np, and so 
Zc{P\y) is expected to take values around (1 ~ pY^e'^"^^^ . The more complicated and more 
interesting question is how does Ze{P\y) behave, and here the treatment will be very similar 
to that of the REM. 

Given y, define ^y{d) as the number of incorrect codewords whose Hamming distance 
from y is exactly d. Thus, 

n 

ZMy) = {l-pf^Y.^y{d)-e-^''. (209) 

d=0 



^®This one is still independent, but it has a different distribution, and hence will be handled separately. 
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Just like in the REM, here too the enumerator ^y{d) is the sum of an exponential number, 
gni?^ of binary i.i.d. RV's: 



%(^)= E X{d{x,y)=d}. 

xec\{Xo} 



(210) 



According to the method of types, the probability of a single 'success' {d{X,y) = n6} is 
given by 



PT{d{X,y) = n6} 



. e 



nh2(5) 



exp{-n[ln2 - h2{6)]}. 



(211) 



So, just like in the REM, we have an exponential number of trials, e"'^, each one with an 
exponentially decaying probability of success, e-"[in2-/i2{5)]_ already know how does this 
experiment behave: It depends which exponent is faster. If i? > ln2 — /;,2(5), we will typically 
see about exp{r;,[i? + h2{5) — In 2]} codewords at distance d = n6 from y. Otherwise, we see 
none. So the critical value of 6 is the solution to the equation 



R + h2i5)-\n2 = 0. 



(212) 



There are two solutions to this equation, which are symmetric about 1/2. The smaller one is 
called the Gilbert-Varshamov (G-V) distancq^ and it will be denoted by 6gv{R) (see Fig. 
[Til) . The other solution is, of course, 6 = 1 — ScviR)- Thus, the condition i? > ln2 — /i2(5) 




ScviR) 0.5 



Figure 14: The Gilbert-Varshamov distance as the smaller solution to the equation R + h2{S) 
ln2 = 0. 
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The G-V distance was originally defined and used in coding theory for the BSC. 
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is equivalent to 5av{R) < S < 1 — 5av{R), ^-nd so, for a typical code in the ensemble: 
Ze{l3\y) « exp{n[i? + /i2(<^)-ln2]}-e-^-^"'^ 

S=Sgv{R) 

l-SaviR) 

= (l-p)/5'^e"(^-^°') ■ exp{n[h2i6)-/3JS]} 

5=5gv{R) 

5gv(^)<'5<1-<5gi'(R) J 

Now, similarly as in the REM, we have to maximize a certain function within a limited 
interval. And again, there are two phases, corresponding to whether the maximizer falls 
at an edgepoint (glassy phase) or at an internal point with zero derivative (paramagnetic 
phase). It is easy to show (exercise: fill in the details) that in the paramagnetic phase, the 
maximum is attained at 

and then 

0(/3) = i?-ln2 + ln[p^ + (l-p)'^]. (214) 
In the glassy phase, 6* — Sgv{R) and then 

m = f3[5Gv{R) Inp + (1 - 5gv{R)) ln(l - p)], (215) 

which is again, linear in (3 and hence corresponds to zero entropy. The boundary between 
the two phases occurs when (3 is such that 5gv{R) — Pi3, which is equivalent to 

^ ^ ^'^^^ ^ ln[(l-p)/p] ■ ^''^^ 

So ^ < /3c{R) is the paramagnetic phase of and /3 > /3c{R) is its glassy phase. 

But now we should remember that Ze is only part of the partition function and it is time 
to put the contribution of Z^. back into the picture. Checking the dominant contribution of 
Z — Zg + Zc as a. function of (3 and R, we can draw a phase diagram, where we find that there 
are actually three phases, two contributed by Zg, as we have already seen (paramagnetic and 
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glassy), plus a third phase - contributed by Zc, namely, the ordered or the ferromagnetic 
phase, where dominates (cf. Fig. [T^ . which means reliable communication, as the correct 
codeword a^o dominates the partition function and hence the posterior distribution. The 
boundaries of the ferromagnetic phase designate phase transitions from reliable to unreliable 
decoding. 



Both the glassy phase and the paramagnetic phase correspond to unreliable communi- 
cation. What is the essential difference between them? As in the REM, the difference is 
that in the glassy phase, Z is dominated by a subexponential number of codewords at the 
'ground-state energy', namely, that minimum seen distance of nScviti), whereas in the para- 
magnetic phase, the dominant contribution comes from an exponential number of codewords 
at distance np^. In the glassy phase, there is seemingly a smaller degree of uncertainty since 
H{X\Y) that is induced from the finite-temperature posterior has zero entropy. But this is 
fictitious since the main support of the posterior belongs to incorrect codewords. This is to 
say that we may have the illusion that we know quite a lot about the transmitted codeword, 
but what we know is wrong! This is like an event of an undetected error. In both glassy 
and paramagnetic phases, above capacity, the ranking of the correct codword, in the list of 
decreasing Pi^{x\y), is about e^^^~'^\ 

Exercise: convince yourself that the phase diagram is as depicted in Fig. [15] and find the 
equations of the boundaries between phases. Note that the triple point is (C, 1) where 



T = 1//3 




C 



Figure 15: Phase diagram of the finite-temperature MAP decoder. 
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C = ln2 — /i2(p) is the channel capacity. Also, the ferro-glassy boundary is the vertical 
straight hne R — C. What does this mean? □ 



5.3 Random Coding Exponents 

It turns out that these findings are relevant to ensemble performance analysis of codes. 
This is because many of the bounds on code performance include summations of P^{y\x) 
(for some which are exactly the partition functions that we work with in the foregoing 
discussion. These considerations can sometimes even help to get tighter bounds. We will now 
demonstrate this point in the context of the analysis of the probability of correct decoding 
above capacity. 

First, we have 



The expression in the square brackets is readily identified with the partition function, and 

we note that the combination oi R > C and /3 ^ oo takes us deep into the glassy phase. 
Taking the ensemble average, we get: 



At this point, the traditional approach would be to insert the expectation into the square 
brackets by applying Jensen's inequality (for /3 > 1), which would give us an upper bound. 

Instead, our previous treatment of random code ensembles as a REM-like model can give 
us a hand on exponentially tight evaluation of the last expression, with Jensen's inequality 




M 




y 





(217) 
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being avoided. Consider the following chain: 







1/13' 














.xec 







-PJd 



.d=0 



1//3' 



max VL'uid)e 

0<d<n ^ 



1//3' 



= {I -PTE 

\0<d<n 

{n 
d=0 

n 

= (i-pr^i;{[%(d)]V/^}-e 

d=0 

= {l-prma^E{[ny{d)Y/^}.e 



Jd 



Jd 



Jd 



-Jd 



Thus, it boils down to the calculation of (non-intcgcr) moments of Qy{d,). At this point, we 
adopt the main ideas of the treatment of the REM, distinguishing between the values of S 
below the G-V distance, and those that are above it. Before we actually assess the moments 
of Qy{d), we take a closer look at the asymptotic behavior of these RV's. This will also 
rigorize our earher discussion on the Gaussian REM. 

For two numbers a and b in [0, 1], let us define the binary divergence as 



D(a\\b) = aln^ + (1 - a) In — y. 

b 1 — 



(218) 



Using the inequality 



ln(l + x) = - In ( 1 - 



X 



1 + x 



> 



X 



1 + x' 
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we get the following lower bound to D{a\\b): 



D(a\\b) = aln^ + (l-a)ln^ 

1 — 

1 \, b — a 

= a In - + (1 - a) In 1 + 



b ^ V 1-^ 

, a , , (6 — a)/(l — 6) 
> aln- + (l-a) ^ ^'^ ' 



b ' ' l + (6-a)/(l-6) 

= aln- + o — a 
b 

> a(ln^-l) 

Now, as mentioned earlier, ^y{d) is the sum of e"^ i.i.d. binary RV's, i.e., Bernoulli RV's 
with parameter g-"[in2-fe2('5)]_ Consider the event ^y{d) > e"'^, A > 0, which means that 
the relative frequency of 'successes' exceeds ^ — e~'^^^~^K Then this is a large deviations 
event if e-"(^-^) > e-"!^'^ 2-/^2 (5)1^ that is, 

A> R + h2{S) -ln2. (219) 

Using the Chernoff bound (exercise: fill in the details), one can easily show that 

PT{Qy{d) > e"^} < exp{-e"-^D(e-"(-^-^)||e-"[''^2-'*^(*)l)}. (220) 

Note: we have emphasized the use of the Chernoff bound as opposed to the method of types 
since the method of types would introduce the factor of the number of type classes, which 
is in this case (e"^ + 1). Now, by applying the above lower bound to the binary divergence, 
we can further upper bound the last expression as 

Pr{ny{d) > e"^} < exp{-e"^ • e""^^"^^ • (n[ln2 - R - h2{5) +A]- 1)} 
= exp{-e"^ ■ (n[ln2 - R - h2{5) + A] - 1)} 

Now, suppose first that 5gv{R) < S < 1 — 5gv{R), and take A — R -\- h2{S) — ln2 + e, 
where e > may not necessarily be small. In this case, the term in the square brackets 
is e, which means that the right-most side decays doubly-exponentially rapidly. Thus, for 
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Sgv{.R) < S < 1 — Sgv{R), the probability that ^y{d) exceeds E{Qy{d)} ■ e"''' decays double- 
exponentially fast with n. One can show in a similar manner (exercise: please doj^ that 
PT{Qy{d) < E{Qy{d)} ■ e~"^} decays in a double exponential rate as well. Finally, consider 
the case where 6 < 6gv{R) or 5 > 1 — 6gv{R), and let A = 0. This is also a large deviations 
event, and hence the above bound continues to be valid. Here, by setting A = 0, we get an 
ordinary exponential decay: 

Pr{%(rf) > 1} < e-n[ln2-i?-/«2(<5)]_ (221) 

Now, after having prepared these results, let's get back to the evaluation of the moments 
of Qy{d). Once again, we separate between the two ranges of 6. For 6 < 6gv{R) or 
6 > 1 — 6gv{R), "we have the following: 

E{[Qy{d)Y/^} = 0^/^ ■ Pr{%(rf) = 0} + e"-°/^ ■ Pr{l < Qy{d) < e"^} + double-exp. terms 
= e"-°/^-Pr{%(c/) > 1} 

Thus, in this range, E{[ny{d)Y^'^} = e-'^[i'^2-i?-/.2(5)] independently of (3. On the other hand 
in the range 6gv{R) < S < 1 — 6gv{R), 

E{[Qy{d)Y/^} = (e"[^+'^2('5)-l-2])l//3 . p^|gn[i?+h2(5)-ln2-.] < ^^(^^ < ^n[R+h,iS)-ln2+e]^ ^ 

+double-exp. terms 

^ ^n[R+h.2{S)-ln2]/l3 

since the probability Pr{e"[-f^+'^2W-in2-e] < ^^^^^-^ < ^niR+h2(S)-in2+e]^ ^^^^^ ^^^^^ double- 
exponentially rapidly. So to summarize, we have shown that the moment of ^y{d) undergoes 
a phase transition, as it behaves as follows: 

jpnn _ / e"[^+'^^W-i°2] s < s^^^r) or 5 > 1 - 5GviR) ....^ 

E{Py[d}\ I - I <6<1- 6gv{R) ^ ' 
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This requires a slighly different lower bound to the binary divergence. 
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Finally, by plugging these moments back into the expression of Pc (exercise: fill in the 
details), and taking the limit ^ — > oo, we eventually get: 



hm E 

B-^oo 



e-""^^ (223) 



.xec 

where Fg is the free energy of the glassy phase, i.e., 

Fg = Sgv{R) In - + (1 - 5gv{R)) In (224) 
p 1-p 

and so, we obtain a very simple relation between the exponent of Pc and the free energy of 
the glassy phase: 



p — p-"^s 



M 

y 

ex.p{n{\n2- R - Fg)} 



= exp{n[ln2 -R + Sgv{R) \np + (1 - 6gv{R)) ln(l - p)]} 
= exp{n[h2{dGv{R)) + 5gv{R) Inp + (1 - 5gv{R)) ln(l - p)]} 

^ ^-nD{SGv{R)\\p) 

The last expression has an intuitive interpretation. It answers the following question: what 
is the probability that the channel would fiip less than n5Gv{R) bits although p > 5gv{RV- 
This is exactly the relevant question for correct decoding in the glassy phase, because in 
that phase, there is a "belt" of codewords "surrounding" y at radius n5Gv{R) - these are 
the codewords that dominate the partition function in the glassy phase and there are no 
codewords closer to y. The event of correct decoding happens if the channel flips less than 
n5Gv{R) bits and then Xq is closer to y more than all belt-codewords. Thus, Xq is decoded 
correctly. 

One can also derive an upper bound on the error probability ai R < C . The partition 
function Z[f5\y) plays a role there too according to Gallager's classical bounds. We will not 
delve now into it, but we only comment that in that case, the calculation is performed in 
the paramagnetic regime rather than the glassy regime that we have seen in the calculation 
of Pc- The basic technique, however, is essentially the same. 
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We will now demonstrate the usefulness of this technique of assessing moments of distance 
enumerators in a certain problem of decoding with an erasure option. Consider the BSC 
with a crossover probability p < 1/2, which is unknown and one employs a universal detector 
that operates according to the following decision rule: Select the message m if 

f,-n^h{XmS)y) 

^ > (225) 

where /3 > is an inverse temperature parameter and h(x (B y) is the binary entropy 
pertaining to the relative number of I's in the vector resulting from bit-by-bit XOR of x 
and y, namely, the binary entropy function computed at the normalized Hamming distance 
between x and y. If no message m satisfies (I225p . then an erasure is declared. 

We have no optimality claims regarding this decision rule, but arguably, it is a reasonable 
decision rule (and hence there is motivation to analyze its performance): It is a universal 
version of the optimum decision rule: 

Decide on m erase otherwise. (226) 

The minimization of h{Xm © y) among all codevectors {a^^}, namely, the minimum con- 
ditional entropy decoder is a well-known universal decoding rule in the ordinary decoding 
regime, without erasures, which in the simple case of the BSC, is equivalent to the maxi- 
mum mutual information (MMI) decoder and to the generalized likelihood ratio test (GLRT) 
decoder, which jointly maximizes the likelihood over both the message and the unknown pa- 
rameter. Here we adapt the minimum conditional entropy decoder to the structure proposed 
by the optimum decoder with erasures, where the (unknown) likelihood of each codeword 
Xm is basically replaced by its maximum Q-'n-h{Xm®y) ^ i^^^ with an additional degree of free- 
dom of scaling the exponent by /3. The parameter /3 controls the relative importance of the 
codeword with the second highest score. For example, when /3 — )■ cxd^J only the first and 
the second highest scores count in the decision, whereas if /3 — )■ 0, the differences between 
the scores of all codewords are washed out. 



^As /3 varies it is plausible to let T scale linearly with (3. 
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To demonstrate the advantage of the proposed analysis technique, we will now apply it 
in comparison to the traditional approach of using Jensen's inequality and supplementing 
an additional parameter p in the bound so as to monitor the loss of tightness due to the use 
of Jensen's inequality. Let us analyze the probability of the event £i that the transmitted 
codeword does not satisfy fl225p . We then have the following chain of inequalities, where 
the first few steps are common to the two analysis methods to be compared: 



Pr{^i} 



1 

m=l y 

M 



< 



M 
1 



M 

m=l y 



M 



y\Xr. 



Q-nl3h{X,„®y) 

nT -nl3h{X^,(By) 
g-nl3h(Xm®y) 

. nl3sh{Xm®y) . 



> 1 



=1 y 



^ ^-nf}hiX^,®y) 
ml 



(227) 



Considering now the ensemble of codewords drawn indepedently by fair coin tossing, we 
have: 



Pr{£i} < e"^^J]^;|p(?/|Xi)-exp[nM(^i©2/)]}-^^ 

y 



^exp[-n/3/i(X^©y)] 

m>l 



A nsT ' 



e"^^5^A(2;)-i?(y) (228) 

y 

The computation of A{y) is as follows: Denoting the Hamming weight of a binary sequence 
z by w(z), we have: 

•Kx®y) 

exp[?T,/3s/i(a; © y)] 



1 — p 

2 

1 — p 
2 

1 — p 



exp 



z 



n ( wiz) In h Psh(z) 

1 — p 



exp 



n f3sh{6)-6\n 



exp 



72 max I (1 + (3s)h(5) — Sin 



1 — p 
p 

1 — p 
p 



(229) 



It is readily seen by ordinary optimization that 

1 — p 



max 

5 



[1 + Ps)h{6)-6\n■ 



p 



1 + (3s) In + (1 _ _ in(i _ p) (230) 



114 



and so upon substituting back into the the bound on Pr{£^i}, we get: 



Pr{^i} < exp [n (sT+ (1 + /3s) In + (i _ _ 1^2)] ■ (231) 

y 

It remains then to assess the exponential order of B[y) and this will now be done in two 
different ways. The first is Forney's way of using Jensen's inequality and introducing the 
additional parameter p, i.e., 

B{y) = eI \ Y,^M-nPkX^®y) 



_m>l 



< 



e\[Y^ eM-n(3sh{Xm ® y)/p] 



< e"''^(^;{exp[-nM(X^©2/)/p]} 



< s/p < 1 
P<1 



(232) 



where in the second line we have used the following inequalitjo for non-negative {aj} and 
9e[0,l]: 

^E^'- (233) 



Now, 



^;|exp[-nM(^m©2;)/p]} = 2-"^exp[-n/3s/i(z)/p] 



,nh{5) . -nl3shiS)/p 



exp[n([l-/3s/p]+-l)ln2]. 



(234) 



where [u]^ = max{M,0}. Thus, we get 



Biy) < exp(n[p(i? - ln2) + [p - /3s]+]), 



(235) 



^^To see why this is true, think of pi = ai/{J2i'''i) ^ probabihties, and then > pi, which imphes 
J2iPi — J2iPi = 1- The idea behind the introduction of the new parameter p is to monitor the possible loss 
of exponential tightness due to the use of Jensen's inequality. If p = 1, there is no loss at all due to Jensen, 
but there is maximum loss in the second line of the chain, li p = s, it is the other way around. Hopefully, 
after optimization over p, the overall loss in tightness is minimized. 
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which when substituted back into the bound on Pr{£^i}, yields an exponential rate of 



E^{R,T) = max {(p-[p-/3s]+)ln2- 



-(1 + ps) In [pV(i+/3^) + (1 _ -pR- sT] 

On the other hand, estimating B{y) by the new method, we have: 



(236) 



B{y) = E 



= E 



^ exp[-n/3/i(X„ e y)] 
^%(n(5)exp[-n/3/i((5)] 



• exp(-n^s/i(5)) 

5 

= ^niR+hiS)-in2] . exp[-n/3s/i(5)] + ^ e"^[^+'*('5)-'"2l . exp[-n/3s/i(5)] 

= + (237) 



where Qr = {S : dcviR) < ^ < ^ — ^Gvi^)}- Now, f/ is dominated by the term 6 = if (3s > 
1 and 6 = doviR) if (3s < 1. It is then easy to see that U = exp[— n(ln2 — i?)(l — [1 — /3s]+)]. 
Similarly, V is dominated by the term 5 = 1/2 if ^<1 and 5 — 5gv{R) if > 1- Thus, 
y = exp[-ns(^[ln2 - R] - R[l - ^]+)]. Therefore, defining 



<j>{R, /3, s) = min{(ln 2 - i?)(l - [1 - /3s]+), s(/3[ln 2-R]-R[l- /3]+)}, 



the resulting exponent is 



(238) 



E^{R, T) = max U{R, /3, s) - (1 + /3s) In W/Ci+z?^) + (i _ _ sT] . (239) 

Numerical comparisons show that while there are many quadruples (p, (3,R,T) for which 
the two exponents coincide, there are also situations where Ei{R,T) exceeds Ei{R,T). To 
demonstrate these situations, consider the values p — 0.1, (3 — 0.5, T — 0.001, and let R 
vary from to 0.06 in steps of 0.01. Table 1 summarizes numerical values of both exponents, 
where the optimizations over p and s were conducted by an exhaustive search with a step 
size of 0.005 in each parameter. In the case of Ei{R,T), where s > is not limited to the 
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R = om 


i? = 0.01 


R = 0.02 


i? = 0.03 


R = 0.04 


R = 0.05 


i? = 0.06 


Ei{R,T) 


0.1390 


0.1290 


0.1190 


0.1090 


0.0990 


0.0890 


0.0790 


Ei{R,T) 


0.2211 


0.2027 


0.1838 


0.1642 


0.1441 


0.1231 


0.1015 



Table 1: Numerical values of Ei{R, T) and Ei{R, T) as functions of i? for p = 0.1, (3 = 0.5, 
and T = 0.001. 



interval [0, 1] (since Jensen's inequality is not used), the numerical search over s was limited 
to the interval [0,5]Ef 

As can be seen (see also Fig. [T6l) . the numerical values of the exponent Ei{R,T) are 
considerably larger than those of Ei{R,T) in this example, which means that the analysis 
technique proposed here, not only simplifies exponential error bounds, but sometimes leads 
also to significantly tighter bounds. 

0.25 1 1 1 1 1 1 1 1 1 1 1 




0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 

R 

Figure 16: Graphs of Ei{R,T) (solid line) and Ei{R,T) (dashed line) as functions of R for 
p = 0.1, T = 0.001 and;g = 0.5. 



There are other examples where these techniques are used in more involved situations. 

It is interesting to note that for some values of R, the optimum value s* of the parameter s was indeed 
larger than 1. For example, at rate R = 0, we have s* = 2 in the above search resolution. 
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and in some of them they yield better performance bounds compared to traditional methods. 
Here is a partial list of papers: 

• R. Etkin, N. Merhav and E. Ordentlich, "Error exponents of optimum decoding for the 
interference channel," IEEE Trans. Inform. Theory, vol. 56, no. 1, pp. 40-56, January 
2010. 

• Y. Kaspi and N. Merhav, "Error exponents of optimum decoding for the degraded 
broadcast channel using moments of type class enumerators," Proc. ISIT 2009, pp. 
2507-2511, Seoul, South Korea, June-July 2009. Full version: available in arXiv:0906.1339. 

• A. Somekh-Baruch and N. Merhav, "Exact random coding exponents for erasure de- 
coding," to appear in Proc. ISIT 2010, June 2010, Austin, Texas, U.S.A. 
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6 Additional Topics (Optional) 

6.1 The REM With a Magnetic Field and Joint Source-Channel 
Coding 

6.1.1 Magnetic Properties of the REM 

Earlier, we studied the REM in the absence of an external magnetic field. The Gaussian 
randomly drawn energies that we discussed were a caricature of the interaction energies in 
the p-spin glass model for an extremely large level of disorder, in the absence of a magnetic 
field. 

We are now going to expand the analysis of the REM so as to incorporate also an external 
magnetic field B. This will turn out to be relevant to a more general communication setting, 
namely, that of joint source-channel coding, where as we shall see, the possible skewedness 
of the probability disitribution of the source (when it is not symmetric) plays a role that is 
analogous to that of a magnetic field. The Hamiltonian in the presence of the magnetic field 
is 

n 

S{s) = -Bj2si + Si{s) (240) 

i=l 

where Si{s) stands for the interaction energy, previously modeled to be A/'(0, ^nJ^) according 
to the REM. Thus, the partition function is now 

s 
s 



, Si 

n 



= E 



m 

A 



^ g-/3£/(S) 
S: m{S)=m 

= 5^Zo(/3,m)-e+"'^^"^ 



m 



where Zo{(3,m) is the partial partition function, defined to be the expression in the square 
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brackets in the second to the last hneo Now, observe that Zo(/3,m) is just hke the par- 
tition function of the REM without magnetic field, except that it has a smaller number 
of configurations - only those with magnetization m, namely, about exp{nh2{{l + m)/2)} 
configurations. Thus, the analysis of Zo{/3,m) is precisely the same as in the REM except 
that every occurrence of the term In 2 should be replaced by /i2((l + m)/2). Accordingly, 



(241) 



with 



max 

\^\<Jy/h2{il+m)/2) 



1 + m 



J 



/3e 



and from the above relation between Z and Zq, we readily have the Legendre relation 

0(/3,5) = max[^(/3,m) + /3mS]. (242) 

m 

For small /3 (high temperature), the maximizing (dominant) m is attained with zero-derivative: 



dm 



ho 



m 



+ 



(5mB 



that is 



which yields 



2 1 + m 



+ = 



(243) 



(244) 



m 



mp{/3, B) = tanh(/35) (245) 

which is exactly the paramagnetic characteristic of magnetization vs. magnetic field (like 
that of i.i.d. spins), hence the name "paramagnetic phase." Thus, plugging m* = tanh(/3i?) 
back into the expression of (p, we get: 

'1 + tanh(/35)\ l3^J^ 



0(/3,fi) = /l2 



4 



+ l3Btanh{l3B). 



(246) 



^^Note that the relation between Zo{(3,m) to Z{l3,B) is similar to the relation between il{E) of the 
microcanonical ensemble to Z{l3) of the canonical one (a Legendre relation in the log domain): we are 
replacing the fixed magnetization m, which is an extensive quantity, by an intensive variable B that controls 
its average. 
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This solution is valid as long as the condition 



/3</3^. 

holds, or equivalently, the condition 



J 



(247) 



^ ^ I l + tanh(/3E) ^ 



4 - V 2 
Now, let us denote by f3c{B) the solution f3 to the equation: 

, /l + tanh(/35)\ 
hi ^ 1. 



(248) 



(249) 



As can be seen from the graphical illustration (Fig. [T7|) . Pc{B) is a decreasing function and 
hence Tc{B) = l/f3c{B) is increasing. Thus, the phase transition temperature is increasing 
with \B\ (see Fig. [T8i) . Below /3 = f3c{B), we are in the glassy phase, where is given by: 




/3c(S2) /3c(Sl)/3c(0) 



Figure 17: Graphical presentation of the solution f3c{B) to the equation \/3 J = /i2((l + 
tanh(/3i?))/2) for various values of B. 



(j){/3, B) = max 



/3 ■ max 



JW/i2 ( — r — 1 +mB 



(250) 



thus, the maximizing m does not depend on /3, only on 5. On the other hand, it should be 
the same solution that we get on the boundary (3 = (3c{B), and so, it must be: 



m* = mg{B) = tanh(S/3,(5)). 



(251) 
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T 




glassy 



B 



Figure 18: Phase diagram in the B-T plane. 



Thus, in summary 




l+mp(l3,B) 
2 



(252) 



l3JJh 




2 



In both phases B ^ imphes m* — )■ 0, therefore the REM does not exhibit spontaneous 
magnetization, only a glass transition, as described. 

Finally, we mention an important parameter in the physics of magnetic materials - the 
weak-field magnetic susceptibility, which is defined as x = ^g-|B=o- If can readily be shown 
that in the REM case 



The graphical illustration of this function is depicted in Fig. [191 The 1/T behavior for high 
temperature is known as Curie's law. As we heat a magnetic material up, it becomes more 
and more difficult to magnetize. The fact that here x an upper limit of 1/Tc(0) follows 
from the random interactions between spins, which make the magnetization more difficult 
too. 

6.1.2 Relation to Joint Source— Channel Coding 

We now relate these derivations to the behavior of joint source-channel coding systems. The 
full details of this part are in: N. Merhav, "The random energy model in a magnetic field 




(253) 
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Tc{0) 




Figure 19: x vs. 


T. 



T 



and joint source-channel coding," Physica A: Statistical Mechanics and Its Applications, vol. 
387, issue 22, pp. 5662-5674, September 15, 2008. 

Consider again our coded communication system with a few slight modifications (cf. Fig. 
I20l) . Rather than e"^ equiprobable messages for channel coding, we are now talking about 
joint source-channel coding where the message probabilities are skewed by the source prob- 
ability distribution, which may not be symmetric. In particular, we consider the following: 
Suppose we have a vector s G { — 1,-|-1}^ emitted from a binary memoryless source with 
symbol probabilities q = Prj^j = +1} = 1 — Prj^j = —1}. The channel is still a BSC with 
crossover p. For every A^-tuple emitted by the source, the channel conveys n channel binary 
symbols, which are the components of a codeword x G {0, 1}*^, such that the ratio 6 = n/N, 
the bandwidth expansion factor, remains fixed. The mapping from s to a; is the encoder. As 
before, we shall concern ourselves with random codes, namely, for every s G {—1, +1}^, we 
randomly select an independent codevector x{s) G {0, 1}" by fair coin tossing, as before. 
Thus, we randomly select 2^ codevectors, each one of length n = NO. As in the case of pure 



s 




x{s) 




y 




s 




encoder 


P{y\x) 


decoder 













Figure 20: Block diagram of joint source-channel communication system. 
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channel coding, we consider the finite-temperature posterior: 
with 

^(^|2/) = E[^(«)^(2/I^(*))]^ (255) 
s 

corresponding to the finite-temperature decoder: 

Si = argmax V [P{s)P{y\x{s))f . (256) 

S: Si=s 

Once again, we separate the contributions of Zc{(3\y) = [P{so)P{y\x{sQ))]'^ , Sq being the 
true source message, and 

ZeW\y) = E [Pis)Piy\^is))r- (257) 

S^So 

As we shall see quite shortly, behaves like the REM in a magnetic field given hj B = 
|ln Accordingly, we will henceforth denote Ze{f3) also by Ze{f3,B), to emphasize the 
analogy to the REM in a magnetic field. 

To sec that Z,.(i3,B) behaves like the REM in a magnetic field, consider the follow- 
ing: first, denote by Ni{s) the number of +l's in s, so that the magnetization, m{s) = 
j^[^iLi^{si — +1} — Sill 1{ Si — pertaining to spin configuration s, is given by 



m(s) = 2A^i(s)/A^ - 1. Equivalently, Ni{s) = A^(l + m(s))/2, and then 

/ X Nil+m{S))/2 



/ N Nm{S))/2 

= [g(l - g)]^/^e^-(*)^ 
where B is defined as above. By the same token, for the binary symmetric channel we have: 
P{y\x) = p<^«(^'2/)(i - pY-dH{x,y) ^y^-Jdu^xy) ^258) 
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where J — In ^ and dnix, y) is the Hamming distance, as defined earher. Thus, 



E 



,-0Hi/p{y\x{s))] 



X{S): m{S)=m 



E 



l3JdH{x(s),y) 



JNmB 



[g(i-g)r/^(i-prE, 

' X{S): m{S)=m 
m 

The resemblance to the REM in a magnetic field is now self-evident. In analogy to the 
above analysis of the REM, Zo[P, m) here behaves like in the REM without a magnetic field, 
namely, it contains exponentially e^K{^+"^)/'^) — gn/i((i+m)/2)/0 ^gpj^g^ with the random energy 
levels of the REM being replaced now by random Hamming distances {dH{x{s),y)} that 
are induced by the random selection of the code {x(s)}. Using the same considerations as 
with the REM in channel coding, we now get (exercise: fill in the details): 



//o N A \nZo{f3,m\y) 
ip[p,m) — lim 

n->oo n 



— max 

<S<1—S„ 



\h2 (^^)+/i2(<^)-ln2-/3J(5 



= ^GV { (^-Y~ 



Ih2{^)+h2{pp)-ln2-pjp^ pp>5,, 

-/3J6m P/3 < S, 



where again. 



P/3 



+ (1 -py 



The condition Pi3 > Sm is equivalent to 

(3 < /3o(m) 

Finally, back to the full partition function: 



A 1 , 1 — 5n 



J 



In 



(259) 



(260) 



</)(/3,S)= lim - In 

n— >-oo iV 



5^Zo(/3,m|2/)e^'3^^ 



= max[^V'(/5, m) + PmB] . (261) 



For small enough (3, the dominant m is the one that maximizes [/i2((l + m)/2) + (3mB], 
which is again the paramagnetic magnetization 



m 



mp{f3,B) = tanh(/35). 



(262) 
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Thus, in high decoding temperatures, the source vectors {s} that dominate the posterior 
Pis{s\y) behave hke a paramagnet under a magentic field defined by the prior i? = | In 
In the glassy regime, similarly as before, we get: 



m 



mg{B) = tanh(S/3e(5)) 



(263) 



where this time, (3c{B), the glassy-paramagnetic boundary, is defined as the solution to the 
equation 



In 2 - h2{pf}) = 



1 /l + tanh(^5)\ 



(264) 



The full details are in the paper. Taking now into account also Zc, we get a phase diagram 
as depicted in Fig. [211 Here, 



5n = ^ In ■ 



2 1 - g* 
where q* is the solution to the equation 

h^iq) = e[ln2 - h^ip)], 

namely, it is the boundary between reliable and unreliable communication. 



(265) 



(266) 
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Figure 21: Phase diagram of joint source-channel communication system. 
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6.2 The Generalized Random Energy Model (GREM) and Hier- 
archical Coding 



In the mid-eighties of the previous century, Derrida extended the REM to the generahzed 
REM (GREM), which has an hierarchical tree sturcture to accommodate possible correlations 
between energy levels of various configurations (and hence is somewhat closer to reality). It 
turns out to have direct relevance to performance analysis of codes with a parallel hierarchical 
structure. Hierarchicial structured codes are frequently encountered in many contexts, e.g., 
tree codes, multi-stage codes for progressive coding and successive refinement, codes for the 
degraded broadcast channel, codes with a binning structure (like in G-P and W-Z coding 
and coding for the wiretap channel), and so on. This part is based on the following papers: 

• B. Derrida, "A generalization of the random energy model which includes correlations 
between energies," J. de Physique - Lettres, vol. 46, L-401-107, May 1985. 

• B. Derrida and E. Gardner, "Solution of the generalised random energy model," J. 
Phys. C: Solid State Phys., vol. 19, pp. 2253-2274, 1986. 

• N. Merhav, "The generalized random energy model and its application to the statistical 
physics of ensembles of hierarchical codes," IEEE Trans. Inform. Theory, vol. 55, no. 
3, pp. 1250-1268, March 2009. 

We begin from the physics of the GREM. For simplicity, we limit ourselves to two stages, 
but the discussion and the results extend to any fixed, finite number of stages. The GREM 
is defined by a few parameters: (i) a number < -Ri < ln2 and R2 = ln2 — (ii) a 
number < ai < 1 and 02 = 1 — Oi. Given these parameters, we now partition the set of 
2" configurations into e"^^ groups, each having e"^^^ configurations]^^! The easiest way to 
describe it is with a tree (see Fig. [22]) . each leaf of which represents one spin configuration. 
Now, for each branch in this tree, we randomly draw an independent random variable, which 



Later, we will see that in the analogy to hierarchical codes, Ri and R2 will have the meaning of coding 
rates at two stages of a two-stage code. 
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will be referred to as an energy component: First, for every branch outgoing from the root, 
we randomly draw ~A/'(0,ainJ^/2), 1 < i < e"^i. Then, for each branch 1 < j < e"^^, 
emanating from node no. i, 1 <i < e""^^, we randomly draw e-jj ~ A/'(0, a2nJ^/2). Finally, 
we define the energy associated with each configuration, or equivalently, each leaf indexed 
by (i, j), as Eij = e, + e,,,-, l<i< e^^S 1 < j < e"«^ 




Figure 22: The GREM with K = 2 stages. 



Obviously, the marginal pdf of each Eij is A/'(0, nJ^/2), just like in the ordinary REM. 
However, unlike in the ordinary REM, here the configurational energies {Eij} are correlated: 
Every two leaves with a common parent node i have an energy component in common and 
hence their total energies are correlated. 

An extension of the GREM to K stages is parametrized by Xl^i -Rf = In 2 and Xl^i '^^ ~ 
1, where one first divides the entirety of 2"- configurations into e"-^^ groups, then each such 
group is subdivided into e"^^ subgroups, and so on. For each branch of generation no. 
£, an independent energy component is drawn according to Af{0, aenJ^/2) and the total 
energy pertaining to each configuration, or a leaf, is the sum of energy components along 
the path from the root to that leaf. An extreme case of the GREM is where K = n, which 
is referred to as the directed polymer on a tree or a directed polymer in a random medium. 
We will say a few words about it later, although it has a different asymptotic regime than 
the GREM, because in the GREM, K is assumed fixed while n grows without bound in the 
thermodynamic limit. 
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Returning back to the case oi K — 2 stages, the analysis of the GREM is conceptually 
a simple extension of that of the REM: First, we ask ourselves what is the typical number 
of branches emanating from the root whose first-generation energy component, e^, is about 
e? The answer is very similar to that of the REM: Since we have e"'^^ independent trials of 
an experiment for which the probability of a single success is exponentially e~^^/("''^"i), then 
for a typical realization: 

|e| > nJ^/oiRi 

exp |n i?i - ^ {:^y } |e| < nJ^aiRi ^'^^^^ 

Next, we ask ourselves what is the typical number VL2{E) of configurations with total en- 
ergy about El Obviously, each such configuration should have a first-generation energy 
component e and second-generation energy component E — e, for some e. Thus, 



Oi(e) 



deOi(e) • exp < n 



Ro- 



02 



E 



nJ 



(268) 



It is important to understand here the following point: Here, we no longer zero-out the 
factor 



exp < n 



R2-- 

02 



E 



nJ 



(269) 



when the expression in the square brackets at the exponent becomes negative, as we did in 
the first stage and in the REM. The reason is simple: Given e, we are conducting fli{e) ■ e^^^ 
indepenent trials of an experiment whose success rate is 

(270) 




Thus, whatever counts is whether the entire integrand has a positive exponent or not. 

Consider next the entropy. The entropy behaves as follows: 

\nn2{E) 



E{E) 



lim 

n—^oo 



n 



j:o(e) Eo(s)>o 
-oo j:o{e) < 



(271) 



where So(£') is the exponential rate of the above integral, which after applying the Laplace 
method, is shown to be: 



ME) 



ma x 

\e\<+nJ\/aiRi 



R 



1 / e \^ 1 



02 



E 



nJ 



(272) 
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How does the function S(-E) behave hke? 

It turns out that to answer this question, we will have to distinguish between two cases: 
(i) Ri/ai < i?2/fl2 and (ii) Ri/ai > -R2/ct20 First, observe that So(-E) is an even function, 
i.e., it depends on E only via \E\, and it is monotonoically non-increasing in \E\. Solving 
the optimization problem pertaining to Sq, we readily find: 



ln2-(^)^ 1^1 <^i 



where Ei = nJ ^ Ri/ai. This is a phase transition due to the fact that the maximizing e 
becomes an edgepoint of its allowed interval. Imagine now that we gradually increase \E\ 
from zero upward. Now the question is what is encountered first: The energy level E, where 

jumps to — oo, or Ei where this phase transition happens? In other words, is E < Ei 
OT E > Ell In the former case, the phase transition at Ei will not be apparent because 

jumps to — oo before, and that's it. In this case, according to the first line of So(-E), 



In 2 — {E/njy vanishes at E = nJV^n2 and we get: 




j:(E) = { ^nj) i-i - - (273) 

exactly like in the ordinary REM. It follows then that in this case, 0(/3) which is the Legendre 
transform of Tj{E) will also be like in the ordinary REM, that is: 

As said, the condition for this is: 




or, equivalently. 



nJy/\n2 = E<Ei=nJ\— (275) 

ai 



— >ln2. (276) 
ai 



^^Accordingly, in coding, this will mean a distinction between two cases of the relative coding rates at the 
two stages. 
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On the other hand, in the opposite case, E > Ei, the phase transition at Ei is apparent, 
and so, there are now two phase transtions: 

f In2-(A)^ \E\<E, 
E(S)=<^ E,<\E\<E (277) 

[ -oo 1^1 > E 

and accordingly (exercise: please show this): 

m = \ PJV^, + R, + ^ A < ^ < ^2 = f ^ (278) 

The first line is a purely paramagnetic phase. In the second line, the first-generation branches 
are glassy (there is a subexponential number of dominant ones) but the second-generation 
is still paramagnetic. In the third line, both generations are glassy, i.e., a subexponen- 
tial number of dominant first-level branches, each followed by a subexponential number 
of second-level ones, thus a total of a subexponential number of dominant configurations 
overall. 

Now, there is a small technical question: what is it that guarantees that Pi < ^2 whenever 
Ri/ai < In 2? We now argue that these two inequalities are, in fact, equivalent. In a paper 
by Cover and Ordentlich (IT Transactions, March 1996), the following inequahty is proved 
for two positive vectors (ai, . . . , a^) and (6i, . . . , 6„): 

min ^ < ^^p^ < max ^. (279) 



Thus, 



Ri Ri + i?2 Ri I \ 

mm — < < max — , (280) 

ie{i,2} ai a\ + 02 ie{i,2} 



but in the middle expression the numerator is i?i + i?2 = In 2 and the denominator is 
01+02 = 1, thus it is exactly In 2. In other words. In 2 is always in between R\ja\ and 
R^ja^- So Ri/ai < In 2 iff Ri/ai < -R2/02, which is the case where (3i < (32- To summarize 
our findings thus far, we have shown that: 
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Case A: Ri/ai < R2/(i2 ~ two phase transitions: 

rin2 + ^ /3</3i 
m={ /3Jv^+i?2 + ^ /3i</3</32 (281) 

Case B: Ri/ai > R^/di ~ one phase transition, hke in the REM: 

m = ['f^r^ (282) 

We now move on to our coding problem, this time it is about source coding with a fidehty 
criterion. For simphcity, we will assume a binary symmetric source (BSS) and the Hamming 
distortion. Consider the following hierarchical structure of a code: Given a block length n, 

we break it into two segments of lengths rii and n2 = n — rii. For the first segment, we 
randomly select (by fair coin tossing) a codebook C = {Xi, 1 < i < e"^^^}. For the second 
segment, we do the following: For each 1 < i < e^^^^, we randomly select (again, by fair 
coin tossing) a codebook d — {xij, 1 < j < e"^^^}. Now, given a source vector x e {0, 1}", 
segmentized as {x',x"), the encoder seeks a pair (i,j), 1 < i < e'^^^^, 1 < J < £"■^^2^ 
such that d{x', Xi) + d{x", Xij) is minimum, and then transmits i using riiRi nats and j 
- using n2i?2 nats, thus a total of {riiRi + n2-R2) nats, which means an average rate of 
R — XRi + (1 — A)i?2 nats per symbol, where X — ni/n. Now, there are a few questions that 
naturally arise: 

• What is the motivation for codes of this structure ? The decoder has a reduced delay. It 
can decode the first rii symbols after having received the first riiRi nats, and does not 
have to wait until the entire transmission of length {riiRi + n2i?2) has been received. 
Extending this idea to K even segments of length n/K, the decoding delay is reduced 
from n to n/K. In the limit of K — n, in which case it is a tree code, the decoder is 
actually delay less. 

• What is the relation to the GREM? The hierarchical structure of the code is that of 
a tree, exactly hke the GREM. The role of the energy components at each branch is 
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now played by the segmental distortions d{x',Xi) and d{x",Xij). The parameters Ri 
and i?2 here are similar to those of the GREM. 

• Given an overall rate R, suppose we have the freedom to choose \, Ri and R2, such 
that R — XRi + (1 — A)i?2, are some choice better than others in some sense? This is 
exactly what we are going to check out.. 

As for the performance criterion, here, we choose to examine performance in terms of the 
characteristic function of the overall distortion, ^^[exp{— s • distortion}]. This is, of course, 
a much more informative figure of merit than the average distortion, because in principle, 
it gives information on the entire probability distribution of the distortion. In particular, 
it generates all the moments of the distortion by taking derivatives at s = 0, and it is 
useful in deriving Chernoff bounds on probabilities of large deviations events concerning the 
distortion. More formally, we make the following definitions: Given a code C (any block 
code, not necessarily of the class we defined), and a source vector x, we define 

A{x) ^mmd{x,x), (283) 

XgC 

and we will be interested in the exponential rate of 

^{s) = E{exp[-sA{X)]}. (284) 
This quantity can be easily related to the "partition function" : 

Z{/3\x)^J2^~^'^'''^^- (285) 
xec 

In particular, 

E{exp[-sA{X)]} = lim E {[Z{s ■ eiX)^/^ . (286) 

Thus, to analyze the characteristic function of the distortion, we have to assess (noninteger) 
moments of the partition function. 
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Let's first see what happens with ordinary random block codes, without any structure. 
This calculation is very similar the one we did before in the context of channel coding: 



E{[z{s-e\x)Y/^} 



E 



E 



x&c 



sed{x,x) 



d=0 



'Sdd 




-sd 



d=0 



where, as we have already shown in the past: 



^n[R+h2{5)-ln2] ^ < ScviR) OT 6 > 1 - ScviR) 
^n[R+h,{5yin2]/e < 5 < 1 - dcviR) 



(287) 



Note that 6gv{R) is exactly the distortion-rate function of the BSS w.r.t. the Hamming 
distortion. By plugging the expression of E{['[l(d)Y^^} back into that of E{[Z{s ■ 9\X)Y^^} 
and carrying out the maximization pertaining to the dominant contribution, we eventually 
(exercise: please show that) obtain: 



-nu{s,R) 



(2^ 



where 



u{s,R) = ln2 — i?— max [h2{S) — s6] 

S<5gv{R) 

s5gv{R) S <Sr 
v{s,R) S>Sr 



with 



and 



sr = In 



1 - ScviR) 
5gv{R) 



v{s,R) =ln2-/? + s-ln(l + e"). 
The function u{s, R) is depicted qualitatively in Fig. 



(289) 



(290) 



(291) 
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u(s, R) 




Figure 23: Qualitative graph of the function u{s, R) as a function of s for fixed R. 



Let's now move on to the hierarchical codes. The analogy with the GREM is fairly clear. 
Given x, there are about fti(6i) = e"i[^i+'*2{5i)-in2] first-segment codewords {xi} in C at 
distance ni6i from the first segment x' of x, provided that i?i + /?.2(5i) — ln2 > and Qi{6i) = 
otherwise. For each such first-segment codeword, there are about e'*^!^^"'"'*^^'^^''"^"^! second- 
segment codewords {xij} at distance 712^2 from the second segment x" of x. Therefore, for 

l-daviRi) 



gni[iii+/i2(5i)-ln2] _ gn2[il2+fe2((5-A(5i)/(l-A))-ln2] 



Si=5gv{Ri) 



— exp 



< n ■ 



max 

Sie[SGv(.Ri),i-SGv(.Ri)] 



R + Xh2{Si) + {l-X)h2 



1- A 



In analogy to the analysis of the GREM, here too, there is a distinction between two cases: 
Ri> R> R2 and R\ < R < i?2- In the first case, the behavior is just like in the REM: 



i? + /i2(5) -ln2 5 e[5Gv{R)A-5Gv{R)] 
—00 elsewhere 



(292) 



and then, of course, = —u{/3, R) behaves exactly hke that of a general random code, in 
spite of the hierarchical structure. In the other case, we have two phase transitions: 



-v{f3,R) p<mi) 

(f>{P, R) = { -XpSoviRi) - (1 - XHP, R2) PiRi) <P< PiR2) 
-P[X5GviRi) + (1 - >^)Sgv{R2)] P > 



(293) 



The last line is the purely glassy phase and this is the relevant phase because of the limit 

6* ^ that we take in order to calculate '^{s). Note that at this phase the slope is X5gv{Ri) + 
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(1 — X)5av{R2) which means that code behaves as if the two segments were coded separately, 
which is worse that 5cviR) due to convexity arguments. Let's see this more concretely 

on the characteristic function: This time, it will prove convenient to define Q{di,d2) as an 
enumerator of codewords whose distance is di at the first segment and d2 ~ on the second 
one. Now, 



E {Z^/\s ■ e)} ^ E 



.di=0d2=0 



cii=0d2=0 



(294) 



Here, we should distinguish between four types of terms depending on whether or not 5i e 
[Sgv{Ri): 1 — ^Gv{Ri)] and whether or not 62 G [Sgv{R2): 1 — ^gv{R2)]- In each one of these 
combinations, the behavior is different (the details are in the paper). The final results are 
as follows: 



For Ri<R2, 



lim 



— lni;exp{-sA(X)} 
n 



\u{s,Ri) + {l- \)u{s,R2) 



(295) 



which means the behavior of two independent, decoupled codes for the two segments, 
which is bad, of course. 



For Ri>R2, 



lim 

n— )-oo 



lni;exp{-sA(X)} 



n 



u{s, R) Vs < So 



(296) 



where Sq is some positive constant. This means that the code behaves like an unstruc- 
tured code (with delay) for all s up to a certain sq and the reduced decoding delay 
is obtained for free. Note that the domain of small s is relevant for moments of the 
distortion. For Ri = R2, Sq is unlimited. 



Thus, the conclusion is that if we must work at different rates, it is better to use the higher 
rate first. 
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Finally, we discuss a related model that we mentioned earlier, which can be thought of 
as an extreme case of the GREM with K = n. This is the directed polymer in a random 
medium (DPRM): Consider a Cayley tree, namely, a full balanced tree with branching ratio 
d and depth n (cf. Fig. [2U where d = 2 and n = 3). Let us index the branches by a pair 
of integers {i,j), where 1 < i < n describes the generation (with i = 1 corresponding to 
the d branches that emanate from the root), and < j < c?* — 1 enumerates the branches 
of the i-th generation, say, from left to right (again, see Fig. |2^ . For each branch {i,j), 
^ ^ j ^ d^ , 1 < i < n, we randomly draw an independent random variable Sij according 
to a fixed probability function q{e) (i.e., a probability mass function in the discrete case, or 
probability density function in the continuous case). As explained earlier, the asymptotic 
regime here is different from that of the GREM: In the GREM we had a fixed number of 
stages K that didn't grow with n and exponentially many branches emanating from each 
internal node. Here, we have K = n and a fixed number d of branches outgoing from each 
note. 




Figure 24: A Cayley tree with branching factor d = 2 and depth n = 3. 



A walk w, from the root of the tree to one of its leaves, is described by a finite sequence 



{ihji)}i=i, where < Ji < d - 1 and dji < j^+i < dji + d - 1, i = 1,2, {n - 1) 



27 



For 



^''In fact, for a given n, the number j„ alone dictates the entire walk. 
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a given realization of the RV's {sij : i = 1, 2, . . . , n, j = 0, 1, . . . , — 1}, we define the 
Hamiltonian associated with w as S{w) = Yli=i^i,jiJ then the partition function as: 

= J]exp{-^£:(w;)}. (297) 
w 

It turns out that this model is exactly solvable (in many ways) and one can show (see e.g., 
E. Buffet, A. Patrick, and J. V. Pule, "Directed polymers on trees: a martingale approach," 
J. Phys. A: Math. Gen., vol. 26, pp. 1823-1834, 1993) that it admits a glassy phase transition: 

m = hm i^^^ = I « ^ almost surely (298) 



where 



A ln\d ■ E 



W) = (299) 



and /3c is the value of ^ that minimizes (f)o{f3). 

In analogy to the hierachical codes inspired by the GREM, consider now an ensemble of 
tree codes for encoding source n-tuples, ), which is defined as follows: Given 

a coding rate R (in nats/source-symbol), which is assumed to be the natural logarithm of 
some positive integer d, and given a probability distribution on the reproduction alphabet, 
Q = {Q{y)y y ^ 3^}) l^t us draw d = independent copies of Y under Q, and denote them 
by Yi, Y2, . . . , Yd. We shall refer to the randomly chosen set, Ci = {Yi, Y2, . . . , Y^}, as our 
'codebook' for the first source symbol, Xi. Next, for each 1 < ji < d, we randomly select 
another such codebook under Q, C2J1 = {^1,1)^1,2, ■ ■ ■ j^i.d}) for the second symbol, X2. 
Then, for each 1 < ji < d and 1 < J2 < d, we again draw under Q yet another codebook 
^3,iij2 = {^ij2,i? ^i,j2,2i • • • ! ^i.j2,d}! fo^ so on. In general, for each t < n, we 

randomly draw d^~^ codebooks under Q, which are indexed by (ji, J2, ■ ■ ■ ,jt-i), 1 < jfc < d, 
l<k<t-l. 

Once the above described random code selection process is complete, the resulting set of 
codebooks {Ci, j^. . 2 < t < n, 1 < jfc < (i, 1 < fc < t — 1} is revealed to both the 
encoder and decoder, and the encoding-decoding system works as follows: 



138 



• Encoding: Given a source n-tuple X", find a vector of indices {j*, ■ ■ ■> Jn) ^^^^ 
minimizes the overall distortion X]"=i P(^t? ^i,- -,it)- Represent each component 
(based on hj R = Ind nats (that is, log2 d bits), thus a total of nR nats. 

• Decoding: At each time t {1 < t < n), after having decoded (j*, . . . , j^*), output the 
reproduction symbol Yj*^ j*. 

In order to analyze the rate-distortion performance of this ensemble of codes, we now 
make the following assumption: 

The random coding distribution Q is such that the distribtion of the RV p{x, Y) is the same 
for all X & X. 

It turns out that this assumption is fulfilled quite often - it is the case whenever the 
random coding distribution together with distortion function exhibit a sufficiently high de- 
gree of symmetry. For example, if Q is the uniform distribution over y and the rows of the 
distortion matrix {p{x,y)} are permutations of each other, which is in turn the case, for 
example, when X = y is a group and p{x, y) = j{x — y) is a difference distortion function 
w.r.t. the group difference operation. Somewhat more generally, this assumption still holds 
when the different rows of the distortion matrix are formed by permutations of each other 
subject to the following rule: p{x, y) can be swapped with p(x, y') provided that q{y') — q{y)- 

For a given x and a given realization of the set of codebooks, define the partition function 
in analogy to that of the DPRM: 



where the summation extends over all possible walks, w — {ji, ■ ■ ■ ,jn), along the Cayley 
tree. Clearly, considering our symmetry assumption, this falls exactly under the umbrella 
of the DPRM, with the distortions {p{xt,yji,...,jt)} playing the role of the branch energies 
{si.j}- Therefore, ^lnZ„(/3) converges almost surely, as n grows without bound, to 4>{(3), 
now defined as 



n 




(300) 



w 



t=i 




(301) 
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where now 



Thus, for every (xi, 2:2, • • •)) the distortion is given by 

1 " A 1 

hm sup - ^ p{xt, Yj*^,„j* ) = hm sup - 



t=i 



mm 



= hm sup hm sup 

< hm sup hm sup 
^= ' — hm inf (f){Pe) 

>-oo 



t=i 



_hlZn(A) 

n/3e 



max 

/3>0 



ln[i;{e- 



{x,Y) 



}] + R 



where: (i) {/3^}^>i is an arbitrary sequence tending to inhnity, (ii) the almost-sure equahty 
in the above mentioned paper, and (iii) the justification of the inequahty at the third hne is 
left as an exercise. The last equation is easily obtained by inverting the function R{D) in 
its parametric representation that we have seen earlier: 



R{D) = — min min < I3D + p(x) In 



.y&y 



(302) 



Thus, the ensemble of tree codes achieves R{D) almost surely. 
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6.3 Phase Transitions of the Rate— Distortion Function 



The material in this part is based on the paper: K. Rose, "A mapping approach to rate- 
distortion computation and analysis," IEEE Trans. Inform. Theory, vol. 40, no. 6, pp. 
1939-1952, November 1994. 

We have seen in one of the earlier meetings that the rate-distortion function of a source 
P = {p{x), X G X} can be expressed as 



R(D) = — min 

l3>0 



f3D + J2pix) In 5^g(i/)e- 



■I3d{x,y) 



(303) 



where Q = {q{y),y G 3^} is the output marginal of the test channel, which is also the one 
that minimizes this expression. We are now going to take a closer look at this function in 
the context of the quadratic distortion function d{x,y) = {x — y)^. As said, the optimum Q 
is the one that minimizes the above expression, or equivalently, the free energy 



X \ y / 



(304) 



and in the continuous case, summations should be replaced by integrals: 

'^d{x,y) 



f{Q) = -15 I dxp(x)ln ( / dyq{y)e~ 

P J —OO \J — OO 



(305) 



Rose suggests to represent the RV F as a function of [/ ~ unif[0, 1], and then, instead of 
optimizing Q, one should optimize the function y{u) in: 



f{y{-)) ^ j dxp(x)ln (^j^ d^i{u)e 



-pd{x,y{u)) 



(306) 



where u(-) is the Lebesgue measure (the uniform measure). A necessary condition for opti- 
mality^j which must hold for almost every u is: 



+ 00 



dxp{x) 



-(Sd{x,y{u)) 



/o d/i(M')e-M-,yK)) 



dd{x,y{u)) 
dy{u) 



0. 



(307) 



^^The details are in the paper, but intuitively, instead of a function y{u) of a continuous variable u, think 
of a vector y whose components are indexed by u, which take on values in some grid of [0, 1]. In other words, 



think of the argument of the logarithmic function as X)i=o ^ 
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Now, let us define the support of y as the set of values that y may possibly take on. Thus, 
this support is a subset of the set of all points {yo — y{uo)} for which: 



/" 

J — ( 



dxp{x) 



-I3d{x,ya) 



dd{x, y{u)) 



dy{u) 



0. 



(308) 



y{u)=yo 



This is because yo must be a point that is obtained as y{u) for some u. Let us define now 
the posterior: 



q[u\x) 



g-/3d(x,2/(u)) 



/q di^{u')e-^d{x,y{u' 



)) 



Then, 



, . ^ . , ^ dd(x,y(u)) ^ 
dxp{x)q{u\x) ■ — ^ . . — = 0. 

dy[u) 



(309) 



(310) 



But p{x)q{u\x) is a joint distribution p{x,u), which can also be thought of as fj,{u)p{x\u) . 
So, if we divide the last equation by n{u), we get, for almost all u: 



dxp(x\u) — ; ; = 0. 

dy{u) 



(311) 



Now, let's see what happens in the case of the quadratic distortion, d{x, y) — {x — y^. Let 
us suppose that the support of Y includes some interval Xq as a subset. For a given u, y{u) 
is nothing other than a number, and so the optimality condition must hold for every y G lo- 
in the case of the quadratic distortion, this optimality criterion means 

f + OO 



with 



or, equivalently. 



/ 

J — t 



\{x) 



dxp{x)\{x){x — y)e 



0, e Xo 



/q^ d/i(M)e-/^'i(^'f(")) /^^ dyq{y)e-P'i^^^y) ' 
d 



+ 00 



dxp{x)\{x) 



dy 



-/3{x-yy 



0. 



(312) 



(313) 



(314) 



Since this must hold for all y & Iq, then all derivatives of the l.h.s. must vanish within Xq, 
i.e.. 



+00 ^ 

dxp{x)X{x 

■oo 



ay" 



0. 



(315) 
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J — c 



Now, considering the Hermitian polynomials 

ii^{z)^e^^'^S^-^^') (316) 

this requirement means 

/•+00 

^xp{x)\{x)En{x - y)e-^(^-2')' = 0. (317) 

In words: \{x)'p{x) is orthogonal to all Hermitian polynomials of order > 1 w.r.t. the weight 
function e~'^^^. Now, as is argued in the paper, since these polynomials are complete in 

L^{e~^^'^), we get 

■p{x)\{x) — const. (318) 

because H^iz) = 1 is the only basis function orthogonal to all Hn{z), n > 1. This yields, 
after normalization: 

pix) = ^ j\y^{u)e-^^^-y^^^^' = dyq{y)e-^(^-y^' ^Q^N (o, ^) . (319) 

The interpretation of the last equation is simple: the marginal of X is given by the convo- 
hition between the marginal of Y and the zero-mean Gaussian distribution with variance 
D — 1/(2/3) (= kT/2 of the equipartition theorem, as we already saw). This means that X 
must be representable as 

X^Y + Z (320) 

where Z ~ jV ^0, and independent of Y . From the Information Theory course we know 
that this is exactly what happens when R{D) coincides with its Gaussian lower bound, a.k.a. 
the Shannon lower hound. Here is a reminder of this: 

R{D) = h{X)- max h{X\Y) 

E(X~Yy<D 

= h(X)- max h{X -Y\Y) 

E{X-Y)2<D 

> h{X) - max h{X - Y) equality ii{X -Y) ±Y 

E{X-YY<D 

= h{X) - max h{Z) Z ^ X -Y 
Ez^<D 

> h{X) - \- ln(27reD) equality if Z ~ A^(0, D) 

- ^slb(^) 
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The conclusion then is that if the support of Y includes an interval (no matter how small) 
then R{D) coincides with R^i^j^iD). This implies that in all those cases that -RgLB^-^) 
not attained, the support of the optimum test channel output distribution must be singular, 
i.e., it cannot contain an interval. It can be, for example, a set of isolated points. 

But we also know that whenever R{D) meets the SLB for some D — Dq, then it must 
also coincide with it for all D < Dq. This follows from the following consideration: If X can 
be represented as F + Z, where Z ~ A/'(0, Dq) is independent of F, then for every D < Dq, 
we can always decompose Z as Zi + Z2, where Zi and Z2 are both zero-mean independent 
Gaussian RV's with variances Dq — D and D, respectively. Thus, 

X = y + z = (r + Zi) + = r' + ^2 (321) 

and we have represented X as a noisy version of Y' with noise variance D. Whenever X 
can be thought of as a mixture of Gaussians, R{D) agrees with its SLB for all D upto the 
variance of the narrowest Gaussian in this mixture. Thus, in these cases: 

\ >^slbP) D>Do ^ > 

It follows then that in all these cases, the optimum output marginal contains intervals for 
all D < Dq and then becomes abruptly singular as D exceeds Dq. Prom the viewpoint 
of statistical mechanics, this looks like a phase transition, then. Consider first an infinite 
temperature, i.e., /3 = 0, which means unlimited distortion. In this case, the optimum output 
marginal puts all its mass on one point: y = E[X), so it is definitely singular. This remains 
true even if we increase f3 to the inverse temperature that corresponds to Dmax, the smallest 
distortion for which R{D) — 0. If we further increase /3, the support of Y begins to change. 
In the next step it can include 2 points, then 3 points, etc. Then, if there is Dq below which 
the SLB is met, then the support of Y abruptly becomes one that contains one interval at 
least. This point is also demonstrated numerically in the paper. 

An interesting topic for research evolves around possible extensions of these results to 
more general distortion measures, other than the quadratic distortion measure. 
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6.4 Capacity of the Sherrington— Kirkpartrick Spin Glass 



This part is based on the paper: O. Shental and I. Kanter, "Shannon capacity of infinite- 
range spin-glasses," Technical Report, Bar Ilan University, 2005. In this work, the authors 
consider the S-K model with independent Gaussian coupling coefficients, and they count the 
number N{n) of meta-stable states in the absence of magnetic field. A meta-stable state 
means that each spin is in its preferred polarization according to the net field that it 'feels', 
i.e.. 

Si = sgn JijSj^ , i = l,...,n. (323) 

They refer to the limit hm„^oo[ln N{ri)]/n as the capacity C of the S-K model. However, they 
take an annealed rather than a quenched average, thus the resulting capacity is somewhat 
optimistic. The reason that this work is brought here is that many of the mathematical tools 
we have been exposed to are used here. The main result in this work is that 

C = ln[2(l-g(t))]-^ (324) 

where 

Q{t) ^ — / du- e-"'/2 (325) 
and t is the solution to the equation 

t = . (326) 

V2^[i - Q{t)] ^ ' 

The authors even address a slighlty more general question: Quite obviously, the metasta- 
bility condition is that for every i there exists Aj > such that 

A^Si = ^ ] Jij^j- (327) 

j 

But they actually answer the following question: Given a constant what is the expected 
number of states for which there is Xi > K for each i such that AjSj = Jij^r For 
K — > — oo, one expects C — > ln2, and for X — > oo, one expects C — > 0. The case of interest 
is exactly in the middle, where K = 
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Moving on to the analysis, we first observe that for each such state, 



d\i6 ^ Juse - XiSi 



thus 



N{n) 



fOO POD 



I'OO fOO / " / \ \ 



(328) 



(329) 



Now, according to the S-K model, {Ju} are n{n — l)/2 i.i.d. zero-mean Gaussian RV's with 
variance J^/n. Thus, 



N{n) 



27rJ2 



/ 



dJexp 



n 



2J2 



S "'^ i=l \ I 



(330) 



The next step is to represent each Dirac as an inverse Fourier transform of an exponent 



+ CXD 



dwe^'' 



J = V-1 



(331) 



which then becomes 

n \n(n-l)/4 



N{n) 



^^^^^^^^^^ djexp j-^ g 4[ . x: • • • / 



dA X 



_ / n Nn(n-i)/4 r p 
\2'kP) J^n(n-i)/2 Jk 



dA X 



da; 



exp 



(27r)- 

We now use the Hubbard-Stratonovich transform: 



) XI "^^^ + 5Z Ji^i'^i^i + ~ ^iSiK \ (332) 

I, i>i i>e i ) 



(333) 



with a — n/ (2 J^) and b — uJiS£ + ouiSf 



■ ■ ■ / / 7^ n n ^M-{^iSE + c.,..)'^V(2n)}. (334) 
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Next observe that the summand doesn't actually depend on s because each Sj is multiplied 
by an integration variable that runs over IR and thus the sign of Sj may be absorbed by this 
integration variable anyhow (exercise: convince yourself). Thus, all 2" contributions are the 
same as that of s = (+1, . . . , +1): 

/oo POO r J 

Jk"^^ y^^^ne-^-^»n^^Pi-(^^ + ^^)''^V(2n)}. (335) 

Now, consider the following identity (exercise: prove it): 



j_2 

2n 



and so for large n, 
J2 



(336) 



J- 



J2 



J2 



2n 



Y,{^i+^if ~ Y E'^' + —Y.^^''^ ~ Y E^' + ^ E 



(337) 



i>e 



.1=1 



thus 



(338) 



We now use again the Hubbard-Stratonovich transform 

7iR27r 



(339) 



and then, after changing variables Aj JXi and JcUj cui (exercise: show that), we get: 



CO, 



(340) 
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which after changing tj ^Jn — > becomes 



N{n) 



1 


n 


/ dte 


• — 


/2^J 


1 


n 


/ dte 

IR 






1 




/ dte 

IR 






1 


n 


/ dte 

IR 






n 

V27f 


/ dt 
Jn 


exp 1 



dA / duje^^^'-^^-^"/^ 

K/X Ju 



2^ / dAe-(*-^)^/2 
'k/x 

J —oo 



(again, the H-S identity) 

t^t + K/J, A ^> -A + i + K/J 



in. 



— exp n • max 



ln(2(l - Q{t)) 



{t + K/Jf 



Laplace integration 



The maximizing t zeroes out the derivative, i.e., it solves the equation 

e-^V^ _ ^ K 
x/2^[l - Q(t)] ~ J 

which for K = 0, gives exactly the asserted result about the capacity. 



(341) 



148 



6.5 Generalized Temperature, de Bruijn's Identity, and Fisher In- 
formation 



Earlier, we defined temperature by 

1 fdS 



This definition corresponds to equilibrium. We now describe a generalized definition that is 
valid also for non-equilibrium situations, and see how it relates to concepts in information 
theory and estimation theory, like the Fisher information. The derivations here follow the 
paper: K. R. Narayanan and A. R. Srinivasa, "On the thermodynamic temperature of a 
general distribution," arXiv:0711.1460v2 [cond-mat.stat-mech], Nov. 10, 2007. 

As we know, when the Hamiltonian is quadratic £{x) — f x^, the Boltzmann distribution 
is Gaussian: 

P(x) = iexp|-/3.|f^x?| (343) 
and by the equipartition theorem: 

We also computed the entropy, which is nothing but the entropy of a Gaussian vector S{P) — 
^ln(^). Consider now another probabihty density function Q{x), which means a non- 
equilibrium probability law if it differs from P, and let's look also at the energy and the 
entropy pertaining to Q: 




a 



(345) 



S{Q) = k-{-lnQ{X))Q = -k J dxQ{x)lnQ{x). (346) 

In order to define a notion of generalized temperature, we have to define some sort of deriva- 
tive of S{Q) w.r.t. E{Q). This definition could make sense if it turns out that the ratio 



149 



between the response of S to perturbations in Q and the response of E to the same perurba- 
tions, is independent of the "direction" of this perturbation, as long as it is "small" in some 
reasonable sense. It turns out the de Bruijn identity helps us here. 

Consider now the perturbation of X by V5Z thus defining the perturbed version of X 
as Xs = -X" + V5Z, where 5 > is small and Z is an arbitrary i.i.d. zero-mean random 
vector, not necessarily Gaussian, whose components all have unit variance. Let Qs denote 
the density of Xs (which is, of course, the convolution between Q and the density of Z, 
scaled by \/S). The proposed generalized definition of temperature is: 



1 SiQs)-SiQ) 
T s^o E{Qs) - EiQY 



(347) 



The denominator is easy since 



E\\X + yTbZf - E\Xf = 2V5EX^Z + n5^n5 



(348) 



and so, E{Qs) — E{Q) — na5/2. In view of the above, our new definition of temperature 
becomes: 

(349) 



1 A 2A; ^.^ h{X + v^Z) - h{X) _ 2k dh{X + VSZ) 

T na 5->o S na 85 



6=0 



First, it is important to understand that the numerator of the middle expression is positive 
(and hence so is T) since 



S{Qs) = kh{X + Vsz) > kh{X + VSZ\Z) = kh{X) = S{Q). 



(350) 



In order to move forward from this point, we will need a piece of background. A well-known 
notion from estimation theory is the Fisher information, which is the basis for the Cramer- 
Rao bound for unbiased parameter estimators: Suppose we have a family of pdf's {Qe{x)} 
where ^ is a continuous valued parameter. The Fisher info is defined as 



J{0) = Ee 



dlnQejX) 

de 



-f 

J — c 



Qe{x) 



(351) 



Consider now the special case where ^ is a translation parameter, i.e., Qo{x) — Q{x — 0), 
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then 



+00 

-00 
+00 

-00 
+00 



dx 



d_ 

m 

d_ 

dx 



Q{x - 6) 
Q{x - 9) 



dQ{x- 
dx 



dQ{x - 0) 



J{Q) with a shght abuse of notation. 



independently of 6. For the vector case, we define the Fisher info matrix, whose elements 
are 



Q(x) 



dxi 



dxi 



Shortly, we will relate T with the trace of this matrix. 

To this end, we will need the following result, which is a variant of the well-known de 
Bruijn identity, first for the scalar case: Let Q be the pdf of a scalar RV X of finite variance. 
Let Z he a, unit variance RV which is symmetric around zero, and let Xs ^ X + VSZ. Then, 



dh{X + VSZ) 



85 



J{Q) 



(353) 



5=0 



The original de Bruijn identity allows only a Gaussian perturbation Z, but it holds for any 
5. Here, on the other hand, we allow an arbitrary density M{z) of Z, but we insist on 5 — > 0. 
The proof of this result is essentially similar to the proof of the original result, which can be 
found, for example, in the book by Cover and Thomas: Consider the characteristic functions: 

f+OO 

(354) 



/-I- 00 
dxe'^Q{x) 
-00 



and 



+00 



dze''M{z). 



(355) 
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Due to the independence 



r+oo 



= ^x{s) ■ ^2 — ^^l^ii^n being the i-th moment of Z 



odd moments vanish due to symmetry 



Applying the inverse Fourier transform, we get: 



and so, 



dQs{x) 



05 



(5=0 



1 d^Qjx) 1 d^Qsjx) 

2 ' dx^ ^2' dx^ ' 



Now, let's look at the entropy: 



+00 



h{Xs) = - dxQs{x) In Qs{x). 



Taking the derivative w.r.t. S, we get: 



dhjXs) 
85 



+00 



dx 



dQs{x) dQs{x) 



85 



+ 



In Qs{x) 



8 

l'-.L"" 85 



/-t-00 p+oo 
dxQs{x) — dx 
00 J —00 

8Qs{x) 



85 



85 



\n.Qs{x) 



h\Qs{x) 



-L 



o5 



and so. 



95 



/+00 
dx 
■00 



(9(55 (x) 



95 



5=0 



2 d"x 



Integrating by parts, we obtain: 



dh{Xs) 
85 



5=0 



1 dQ{x) 
'2 ' da; 



+00 



Ing(x) 



+ 



-j^ r+co 
2 7-00 



dx 



"■9g(x) 

8x 



(356) 



(357) 



(358) 



(359) 



"+00 -■ A^Qir^ 
lng(x) = -/ dx-- ' •InQ(x). (360) 



(361) 
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The first term can be shown to vanish (see paper and/or C&T) and the second term is 
exactly J{Q)/2. This completes the proof of the (modified) de Bruijn identity. 

Exercise: Extend this to the vector case, showing that for a vector Z with i.i.d. components, 
all symmetric around the origin: 

Putting all this together, we end up with the following generalized definition of temper- 
ature: 

i = A.tr{J(Q)}. (363) 
1 na 

In the 'stationary' case, where Q is symmetric w.r.t. all components of a;, {Ju} are all the 
same quantity, call it J((5), and then 

1 k 



rr. J{Q) (364) 

1 a 



or, equivalently. 



where CRB is the Cramer-Rao bound. High temperature means a lot of noise and this in 
turn means that it is hard to estimate the mean of X. In the Boltzmann case, J{Q) — 
l/Var{X} — a/3 — a/{kT) and we are back to the ordinary definition of temperature. 

Another way to look at this result is as an extension of the equipartition theorem: As we 
recall, in the ordinary case of a quadratic Hamiltonian and in equilibrium, we have: 

{S{X)) = = ^ (366) 

or 

^ = ^. (367) 

In the passage to the more general case, cr^ should be replaced by 1/J{Q) = CRB. Thus, 
the induced generalized equipartition function, doesn't talk about average energy but about 
the CRB: 

- • CRB = — . (368) 
2 2 ^ ' 
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Now, the CRB is a lower bound to the estimation error which, in this case, is a transaltion 
parameter. For example, let x denote the location of a mass m tied to a spring of strength 
mwl and equilibrium location 6. Then, 

S{x)^'^{x-9f. (369) 

In this case, a — mujQ, and we get: 

estimation error energy = — ^ • E{9{X) - Of > — (370) 

where 0{X) is any unbiased estimator of 9 based on a measurement of X. This is to say 
that the generalized equipartition theorem talks about the estimation error energy in the 
general case. Again, in the Gaussian case, the best estimator is 6{x) — x and we are back 
to ordinary energy and the ordinary equipartition theorem. 
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6.6 The Gibbs Inequality and the Log— Sum Inequahty 

In one of our earlier meetings, we have seen the Gibbs' inequahty, its physical significance, 
and related it to the second law and the DPT. We now wish to take another look at the 
Gibbs' inequality, from a completely different perspective, namely, as a tool for generating 
useful bounds on the free energy, in situations where the exact calculation is difficult (see 
Kardar's book, p. 145). As we show in this part, this inequality is nothing else than the log- 
sum inequality, which is used in Information Theory, mostly for proving certain qualitative 
properties of information measures, like the data processing theorem of the divergence, etc. 
But this equivalence now suggests that the log-sum inequality can perhaps be used in a 
similar way that it is used in physics, and then it could perhaps yields useful bounds on 
certain information measures. We try to demonstrate this point here. 

Suppose we have an Hamiltonian S{x) for which we wish to know the partition function 

Z(^)=^e-^^(^) (371) 

X 

but it is hard, if not impossible, to calculate in closed-form. Suppose further that for another, 
somewhat different Hamiltonian, 6o{x), it is rather easy to make calculations. The Gibbs' 
inequality can be presented as a lower bound on In Z{(3) in terms of B-G statistics pertaining 
to Sq. 



In 




> In 


Y^-pSo(X) 


+ f3{So{X)-£{X))„ 


(372) 




. X 




. X 







The idea now is that we can obtain pretty good bounds thanks to the fact that we may 
have some freedom in the choice of So- For example, one can define a parametric family of 
functions £q and maximize the r.h.s. w.r.t. the parameter(s) of this family, thus obtaining 
the tightest lower bound within the family. We next demonstrate this with an example: 

Example - Non-harmonic oscillator. Consider the potential function 

V{z) = Az"^ (373) 
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and so 



Six) = 1^ + Az\ (374) 
2m 



where we approximate the second term by 

'^"f^) = { +CO N > I 
where L is a parameter to be optimized. Thus, 

-1 r+oo r+oo 

Zo ^ ^ dp d^e-^[^°(^)+^'/(2-)l 

J — oo J — oo 

i / dp . e-^^''/^^-) / dz 

i-oo i-L/2 



V 27rmkT 



h 

and so, by the Gibbs inequahty: 



InZ > lnZo + /3{So{X)-S{X))o 



1 1 



> InZo- ■ ^ / d^->lz^ 

-L/2 



kT L 

> In 



-L/2 



Ly/27rmkT 



h 



80kT 



To maximize f{L) we equate its derivative to zero: 

Plugging this back into the Gibbs lower bound and comparing to the exact value of Z 
(which is still computable in this example), we find that 2'approx ~ 0.9 Inexact' which is 
not that bad considering the fact that the infinite potential well seems to be quite a poor 
approximation to the fourth order power law potential V{z) — Az^. 

As somewhat better approximation is the harmonic one: 

Vo{z) ^!I^.z' (377) 
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where now Uq is the free parameter to be optimized. This gives 

"+00 P+CO 



hj_^ J_oo ruvo 27r 

and this time, we get: 



kT \ 


1 






kT\ 


1 




+ 2- 



= In 

A 



2,. ,4 







Maximizing /: 



^ d/ 1 12AkT , (12AA;T)V4 
= = + =^ uj*o = ^ 379 

This time, we get ^approx ~ 0-95.^gxact' ^^^^ approximation is even better. □ 

So much for physics. Let's look now at the Gibbs inequahty shghtly differently. What we 
actually did, in a nutshell, and in different notation, is the following: Consider the function: 

n n 

Z{X) = al'^b^ = Yl ««e"^'"^"'^'''\ (380) 
i=i i=i 

where {aj} and {6,} are positive reals. Since InZ(A) is convex (as before), we have: 

In (^6,) = InZ(l) 



dA 



In ^ aJ + 



A=0 

En 
i=l 



which is nothing but the log-sum inequality, which in IT, is more customarily written as: 

X:aan^> X^a. -In^r^. (381) 

i=l \i=l J 2^i=l"i 

Returning to the form: 
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the idea now is, once again, to lower bound an expression ln(^"^^ bi) which may be hard 
to calculate, by the expression on the l.h.s. which is hopefully easier, and allows a degree 
of freedom concerning the choice of {oj}, at least in accordance to some structure, and 
depending on a limited set of parameters. 

Consider, for example, a hidden Markov model (HMM), which is the output of a DMC 
W{y\x) = Y[t=i^(yt\^t) fed by a first-order Markov process X, governed by Q{x) — 
YYt=iQ{^t\xt-i)- The entropy rate of the hidden Markov process {Yf} does not admit a 
closed-form expression, so we would like to have at least good bounds. Here, we propose an 
upper bound that stems from the Gibbs inequality, or the log-sum inequality. 

The probability distribution of y is 

n 

P(y) = En[^(^*l^*)^(^^l^*-i)]- (383) 
X t=i 

This summation does not lend itself to a nice closed-form expression, but if the t-th factor 
depended only on t (and not also on t — 1) life would have been easy and simple as the sum 
of products would have boiled down to a product of sums. So this motivates the following 
use of the log-sum inequality: For a given y, let's think of x as the index i of the log-sum 
inequality and then 

n 

b{x) = l[[W{yt\xt)Q{xt\xt-i)]. (384) 
t=i 

Let us now define 

n 

a{x) = l[Po{xt,yt), (385) 
t=i 

where Pq is an arbitrary joint distribution over A" x 3^, to be optimized eventually. Thus, 
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applying the log-sum inequality, we get: 
\nP{y) = \n(^b{x)^ 



Ea;«(^) Hb(x)/a{x)] 



1^ En^oKi/,) + 



X t=l 



^ Exnl^PoM ■ ^'''^ 

Now, let us denote Po{y) = 'Ylix^x Poi^fV)' which is the marginal of y under Pq. Then, the 
first term is simply Ylt=i^^-Po{yt)- As for the second term, we have: 



EE 

t=l X 



J2xIlt=iPo(xt,yt) 
Ylt=i Poi^t: Vt) HQ{xt\xt^^)W{yt\xt)/Po{xt, yQ] 



^ 2^ TT" P r ^ " 2^ i^o(a:;t-i,Z/t-i)^'o(a;t,?/t) -In 



Po(a;t-i,|/t-i)Po(a:t,|/t) 



Po{xt,yt) 



= E E 

t=l Xt-l,Xt 



•In 



Poli/i-Oi'od/t) 
= 5Z XI -fo(2^t-i|?/t-i)-Po(a;t|yt) • In 



t=l Xt-l,Xt 



Q{Xt\X,_,)W{yt\Xt) 
Po{Xt,yt) 



Q{xt\xt-i)W{yt\xt) 
Poixt^yt) 

'Q{xt\xt-i)W{yt\xt) 

Po{Xt:yt) 

Yt-i = yt-i,Yt = yt 



where Eq denotes expectation w.r.t. the product measure of Pq. Adding now the first term 
of the r.h.s. of the log-sum inequality, EtLi lnPo(|/t), we end up with the lower bound: 

■Q(X,|AVi)ir(y,jXi)" 



\TiP{y)>Y,Eo{\^ 



Po{Xt\yt) 



Yt-i = yt-i, Yt = yt} = Y1 Hyt-i,yt; Pq)- 

(387) 

At this stage, we can perform the optimization over Pq for each y individually, and then derive 
the bound on the expectation of In P{y) to get a bound on the entropy. Note, however, that 
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y^; Po) depends on y only via its Markov statistics, i.e., the relative frequencies 
of transitions y =^ y' for all y^y' e y. Thus, the optimum Pq depends on y also via 
these statistics. Now, the expectation of A(|/(_i, y^; Pq) is going to be dominated by the 
typical {y} for which these transition counts converge to the respective joint probabilities of 
{Yi_i — y^ Yf — y}. So, it is expected that for large n, nothing will essentially be lost if we 
first take the expectation over both sides of the log-sum inequality and only then optimize 
over Pq. This would give, assuming stationarity: 

//(y") < ■ max^;{A(yo, ^i; ^o)}- (388) 

where the expectation on the r.h.s. is now under the rea/ joint distribution of two consecutive 
samples of {Yn}, i.e., 

P{yo,yi) = XI '^{xo)Q{xi\xo)P{yo\xo)P{yi\xi), (389) 

X0,X1 

where 7r(-) is the stationary distribution of the underlying Markov process {xt}. 
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6.7 Dynamics, Evolution of Info Measures, and Simulation 



The material here is taken mainly from the books by Reif, Kittel, and F. P. Kelly, Reversibility 
and Stochastic Networks, (Chaps 1-3), J. Wiley & Sons, 1979. 

6.7.1 Mcirkovian Dynamics, Global Balance and Detailed Baleince 

So far we discussed only physical systems in equilibrium. For these systems, the Boltzmann- 
Gibbs distribution is nothing but the stationary distribution of the microstate x at every 
given time instant t. However, this is merely one part of the picture. What is missing is the 
temporal probabilistic behavior, or in other words, the laws that underly the evolution of 
the microstate with time. These are dictated by dynamical properties of the system, which 
constitute the underlying physical laws in the microscopic level. It is customary then to 
model the microstate at time i as a random process {Xt}, where t may denote either discrete 
time or continuous time, and among the various models, one of the most common ones is 
the Markov model. In this section, we discuss a few of the properties of these processes as 
well as the evolution of information measures, like entropy, divergence (and more) associated 
with them. 

We begin with an isolated system in continuous time, which is not necessarily assumed 
to have reached (yet) equilibrium. Let us suppose that Xt, the microstate at time t, can 
take on values in a discrete set X. For r, s e A", let 

w r Pr{Xt+5 = s\Xt = r} 

Wrs — hm — ^ r s (390) 

in other words, Pr{Xt^s = s\Xt = r} = Wrs ■ ^ + o{6). Letting Pr{t) = Pr{Xt = r}, it is 
easy to see that 

Pr{t + dt) = J2 Ps{t)Wsrdt + Pr{t) ['^-^2 ^^^^^ ) ' (^^l) 
s^r \ s^r / 

where the first sum describes the probabilities of all possibile transitions from other states to 
state r and the second term describes the probability of not leaving state r. Subtracting Pr{t) 
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from both sides and dividing by dt, we immediately obtain the following set of differential 
equations: 

dP (t) 

= J2[Psit)Wsr - Prit)Wrs], T E X, (392) 

s 

where Wrr is defined in an arbitrary manner, e.g., Wrr = for all r. These equations are 
called the master equationsc^ When the process reaches stationarity, i.e., for all r G A", 
Pr(t) converge to some Pr that is time-invariant, then 

J^l^sWsr - PrWrs] = 0, V r G A". (393) 

s 

This is called global balance or steady state. When the system is isolated (microcanonical en- 
semble), the steady-state distribution must be uniform, i.e., Pr = l/l'^l for all r E X. From 
quantum mechanical considerations, as well as considerations pertaining to time reversibility 
in the microscopic leveljf^ it is customary to assume Wrs = Wsr for all pairs {r, s}. We then 
observe that, not only, J^sl^^Wgr — PrWrs] = 0, but moreover, each individual term in the 
sum vanishes, as 

PsWsr - PrWrs = ^(^^.r " Wrs) = 0. (394) 

This property is called detailed balance, which is stronger than global balance, and it means 
equilibrium, which is stronger than steady state. While both steady-state and equilibrium 
refer to a situation of time-invariant state probabilities {Pr}, a steady-state still allows cyclic 
flows of probability. For example, a Markov process with cyclic deterministic transitions 
1—7-2— 7-3— 7-1— 7-2— 7-3— T--- - is in steady state provided that the probability distribution 
of the initial state is uniform (1/3, 1/3, 1/3), however, the cyclic flow among the states is in 
one direction. On the other hand, in detailed balance {Wrs = Wsr for an isolated system), 
which is equilibrium, there is no net flow in any cycle of states. All the net cyclic probability 
fluxes vanish, and therefore, time reversal would not change the probability law, that is. 



^^Note that the master equations apply in discrete time too, provided that the derivative at the Lh.s. is 
replaced by a simple difference, Pr{t + 1) — Pr{t), and {Wrs} designate one-step state transition probabilities. 

'^'^Think, for example, of an isolated system of moving particles, obeying the differential equations 
md^ri(t)/dt^ = -^('"ilO ^ '''«(^))) * — Ij 2, . . . , n, which remain valid if the time variable t is replaced 

by -t since d^r,{t)/dt'^ = d^ri{-t)/d{-t)'^. 
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{X_t} has the same probabihty law as {Xt}. For example, if {Yf} is a BernouUi process, 
taking values equiprobably in {—1, +1}, then Xf defined recursively by 

Xt+i = {Xt + Yt)modK, (395) 

has a symmetric state-transition probability matrix W, a uniform stationary state distrib- 
tuion, and it satisfies detailed balance. 

6.7.2 Evolution of Information Measures 

Returning to the case where the process {Xt} pertaining to our isolated system has not 
necessarily reached equilibrium, let us take a look at the entropy of the state 

H{Xt) ^ - J2 Pr{t) log Prit). (396) 

r 

We argue that H{Xt) is monotonically non-decreasing, which is in agreement with the second 
law (a.k.a. the H-Theorem). To this end, we next show that 



dH{Xt] 



> 0, (397) 



dt 

where for convenience, we denote dPr{t)/dt by Pr{t). 

= -Y.[PrmOgPr{t)+Pr{t)] 

r 

= -J2Pr{t)^OgPr{t) ^P,(i)=0 

r r 

= - Z E ^sr[Ps{t) - Pr{t)] logPrit)) W,r = W^. 

r s 

r,s 

lj2^sr[Pr{t) - Ps{t)]logP,{t) 



2 

s,r 

1 



= 2 ^ ^-[^^(t) - Pr{t)] ■ [logPsit) - logP.(t)] 



r.s 



> 0. (398) 
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where the last inequahty is due to the increasing monotonicity of the logarithmic function: 
the product [Ps{t) — Prit)] ■ [logP5(^) — log -Pr(^)] cannot be negative for any pair (r, s), as 
the two factors of this product arc either both negative, both zero, or both positive. Thus, 
H{Xt) cannot decrease with time. 

This result has a discrete-time analogue: If a finite-state Markov process has a symmetric 
transition probability matrix, and so, the stationary state distribution is uniform, then H{Xt) 
is a monotonically non-decreasing sequence. 

A considerably more general result is the following: If {Xt} is a Markov process with a 
given state transition probability matrix W — {Wrs} (not necessarily symmetric) and {Pr} 
is a stationary state distribution, then the function 

U{t)^J2Pr-v[^^ (399) 

is monotonically strictly increasing provided that V{-) is strictly concave. To see why this 
is true, we use the fact that Pg — ^j.^rWrs and define Wgr — PrWrs/Ps- Obviously, 
Ylir — 1 for all s, and so, 

Pr{t + 1) _ ^ Ps{t)Wsr _ J2 ^r.Psjt) ^^^^^ 

Pr Pr Ps 

s s 

and so, by the concavity of V{-): 



EE^»»-i'(^) 

r Q \ 5 / 



=t/(t). (401) 



Here we required nothing except the existence of a stationary distribution. Of course in the 
above derivation t + 1 can be replaced by t + r for any positive real r with the appropriate 
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transition probabilities, so the monotonicity of U{t) applies to continuous-time Markov 
processes as well. 

Now, a few interesting choices of the function V may be considered: 

• For V{x) = —xlnx, we have U{t) = —D{P{t)\\P). This means that the divergence 
between {Pr{t)} and the steady state distribution {Pr} is monotonically strictly de- 
creasing, whose physical interpretation could be the decrease of the free energy, since 
we have already seen that the free energy is the physical counterpart of the divergence. 
This is a more general rule, that governs not only isolated systems, but any Markov 
process with a stationary limiting distribution (e.g., any Markov process whose disti- 
bution converges to that of the Boltzmann-Gibbs distribution). Having said that, if 
we now particularize this result to the case where {Pr} is the uniform distribution (as 
in an isolated system), then 

D{P{t)\\P) = log\X\-H{X,), (402) 

which means that the decrease of divergence is equivalent to the increase in entropy, 
as before. The difference, however, is that here it is more general as we only required 
a uniform steady-state distribution, not necessarily detailed balanceHl 

• Another interesting choice of V is V{x) = Inx, which gives U (t) = —D{P\\P{t)). Thus, 

D{P\\P(t)) is also monotonically decreasing. In fact, both this and the monotonicity 

result of the previous item, are in turn, special cases of a more general result concerning 

the divergence (see also the book by Cover and Thomas, Section 4.4). Let {Pr{t)} and 

{P^{t)} be two time-varying state-distributions pertaining to the same Markov chain, 

but induced by two different initial state distributions, {Pr{0)} and {P^'(O)}. Then 

'^"'^For the uniform distribution to be a stationary distribution, it is sufficient (and necessary) that W would 
be a doubly stochastic matrix, namely, J^r^rs = J2r^sr — 1- This condition is, of course, weaker than 
detailed balance, which means that W is moreover symmetric. 
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D[P{t)\\P'{t)) is monotonically non-increasing. This happens because 



Dipmp'it)) = $^p.(t)iog 



> D{P{t + T)\\P'{t + T)) 




(403) 



where the last inequahty follows from the data processing theorem of the divergence: 
the divergence between two joint distributions of {Xt,Xt^T-) is never smaller than the 
divergence between corresponding marginal distributions of Xt^^-. 

• Yet another choice is V{x) = x'"^, where s G [0, 1] is a parameter. This would yield 
the increasing monotonicity of P^^^P^{t), a metric that plays a role in the theory 
of asymptotic exponents of error probabilities pertaining to the optimum likelihood 
ratio test between two probability distributions. In particular, the choice s = 1/2 
yields balance between the two kinds of error and it is intimately related to the Bhat- 
tacharyya distance. Thus, we obtained some sorts of generalizations of the second law 
to information measures other than entropy. 

For a general Markov process, whose steady state-distribution is not necessarily uniform, 
the condition of detailed balance, which means time-reversibility, reads 



both in discrete time and continuous time (with the corresponding meaning of {VTrs})- The 
physical interpretation is that now our system is (a small) part of a large isolated system, 
which obeys detailed balance w.r.t. the uniform equilibrium distribution, as before. A well 
known example of a process that obeys detailed balance in its more general form is an M/M/1 
queue with an arrival rate A and service rate /x (A < /i). Here, since all states are arranged 
along a line, with bidirectional transitions between neighboring states only (see Fig. 



PsW, 



sr 



PrWr 



(404) 
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there cannot be any cyclic probability flux. The steady-state distibution is well-known to 
be geometric 

which indeed satisfies the detailed balance PrX — -Pr+iA* for all r. Thus, the Markov process 
{Xf}, designating the number of customers in the queue at time t, is time-reversible. 

It is interesting to point out that in order to check for the detailed balance property, 
one does not necessarily have to know the equilibrium distribution {Pr} as above. Applying 
detailed balance to any k pairs of states in a cycle, (si, S2), (^2, S3), ■ ■ ■ , {sk, si), and mul- 
tiplying the respective detailed balance equations, the steady state probabilities cancel out 
and one easily obtains 



W W 

^' S1S2 ^' S1SZ 



w w 



w w 



w w 



(406) 



so this is clearly a necessary condition for detailed balance. One can show conversely, that 
if this equation applies to any finite cycle of states, then the chain satisfies detailed balance, 
and so this is also a sufficient condition. This is true both in discrete time and continuous 
time, with the corresponding meanings of {V^rs} (see Kelly's book, pp. 22-23). 

A A A A 




II II II II 

Figure 25: State transition diagram of an M/M/1 queue. 



In the case of detailed balance, there is another interpretation of the approach to equi- 
librium and the growth of U{t). We can write the master equations as follows: 

dPr{t) fPs{t) PM 



dt 



(407) 



where Rsr — {PsWgr) ^ — (PrWrs) ^. Imagine now an electrical circuit where the indices {r} 
designate the nodes. Nodes r and s are connected by a wire with resistance Rsr and every 
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node r is grounded via a capacitor with capacitance Pr (see Fig. [26]). If Pr(t) is the charge 
at node r at time t, then the master equations are the Kirchoff equations of the currents 
at each node in the circuit. Thus, the way in which probabihty spreads across the circuit 
is analogous to the way charge spreads across the circuit and probabihty fluxes are now 
analogous to electrical currents. If we now choose V{x) = —^x"^, then —U(t) = ^ J2r ^^p^i 
which means that the energy stored in the capacitors dissipates as heat in the wires until 
the system reaches equilibrium, where all nodes have the same potential, Pr{t)/Pr = 1, and 
hence detailed balance corresponds to the situation where all individual currents vanish (not 
only their algebraic sum). 





Figure 26: State transition diagram of a Markov chain (left part) and the electric circuit that 
emulates the dynamics of {Pr{t)} (right part). 



We have seen, in the above examples, that various choices of the function V yield various 
'metrics' between {Pr{t)} and {Pr}, which are both marginal distributions of a single symbol. 
What about joint distributions of two or more symbols? Consider, for example, the function 



(40^ 



where V is concave as before. Here, by the same token, J{t) is a 'metric' between the 
joint probability distribution {P{Xq = r, Xf = s)} and the product of marginals {P{Xq = 
r)P{Xt = s)}, namely, it a measure of the amount of statistical dependence between Xq and 
Xt. For V{x) = Inx, we have, of course, J{t) = —I{XQ;Xt). Now, using a similar chain of 
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inequalities as before, we get the non-decreasing monotonicity of J{t) as follows: 
J(f) - VPfX -r X -s X -v) y( PiXo^r)P{X,^s) PjX.^^ ^ u\X, ^ s) \ 

= 5^P(Xo = r, = ^P(Xi = s\Xo = r, = u) x 

r,u s 

f P{X, = r)P{X, = .s. X,^r = u) \ 
V P{Xo ^r, Xt^ s, Xt+, ^u) ) 

< J]P(Xo = r, Xt+,^u) X 

r,u 

This time, we assumed nothing beyond Markovity (not even homogeneity). This is exactly 
the generalized data processing theorem of Ziv and Zakai (J. Ziv and M. Zakai, "On function- 
als satisfying a data-processing theorem," IEEE Trans. Inform. Theory, vol. IT-19, no. 3, 
pp. 275-283, May 1973), which yields the ordinary data processing theorem (of the mutual 
information) as a special case. Thus, we see that the second law of thermodynamics is (at 
least indirectly) related to the data processing theorem via the fact that they both stem from 
some more general principle concerning monotonic evolution of 'metrics' between probability 
distributions defined using convex functions. In a very similar manner, one can easily show 
that the generalized conditional entropy 

E P(X, = r,X, = s).V (^^^^_i^^) (410) 

7*, 5 

is monotonically non-decreasing with t for any concave V . 
6.7.3 Monte Carlo Simulation 

Returning to the realm of Markov processes with the detailed balance property, suppose we 
want to simulate a physical system, namely, to sample from the Boltzmann-Gibbs distribu- 
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tion 

Pr = (411) 

In other words, we wish to generate a discrete-time Markov process {^t}, possessing the 
detailed balance property, whose marginal converges to the Boltzmann-Gibbs distribution. 
This approach is called dynamic Monte Carlo or Markov chain Monte Carlo (MCMC). How 
should we select the state transition probability matrix W to this end? Substituting Pr — 
e~^^^ / Z[I3) into the detailed balance equation, we readily see that a necessary condition is 

!!z£_g-/3(fi.-£.)_ (412) 

The Metropolis algorithm is one popular way to implement such a Markov process in a rather 
efficient manner. It is based on the concept of factoring Wrs as a product Wrs — Crs-^rsi 
where Crs is the conditional probability of selecting X^+i = s as a candidate for the next state, 
and designates the probability of acceptance. In other words, we first choose a candidate 
according to C, and then make a final decision whether we accept this candidate or stay in 
state r. The Metropolis algorithm pics C to implement a uniform distribution among n states 
'close' to r (e.g., ffipping one spin of a n-spin configuration). Thus, Wrs/Wgr — Ars/Agr, 
and so, it remains to choose A such that 

The Metropohs algorithm defines 

Ars^l , A (414) 

[ 1 otherwise ^ ^ 

In simple words, the algorithm works as follows: Given that Xt = r, first randomly select 

one candidate s for Xt+i among n possible (neighboring) states. If Eg < Er always accept 

Xt+i — s as the next state. If Eg > E^., then randomly draw a RV Y e Unif[0, 1]. If 

Y < e~P{Es-Er)^ then again, accept Xt+i — s as the next state. Otherwise, stay in state r, 

i.e., Xt^i — r. To see why this choice of A works, observe that 

^sr L e-l^('^r-Es) ^ 

There are a few nice things about this algorithm: 
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• Energy differences between neighboring states, Eg — Er, are normally easy to calculate. 
If r and s differ by a single component of the microstate x, and the if the Hamiltonian 

structure consists of short-range interactions only, then most terms of the Hamiltonian 
are the same for r and s, and only a local calculation is required for evaluating the 
energy difference. 

• Calculation of Z{I3) is not required, and 

• Chances are that you don't get stuck in the same state for too long. 

The drawback, however, is that aperiodicity is not guaranteed. This depends on the Hamil- 
tonian. 

The heat bath algorithm (a.k.a. Glauber dynamics) alleviates this shortcoming and al- 
though somewhat slower than Metropolis to equilibrate, it guarantees all the good properties 
of a Markov chain: irreducibility, aperiodicity, and convergence to stationarity. The only 
difference is that instead of the above choice of Ars, it is redefined as 



which is also easily shown to satisfy the detailed balance condition. The heat bath algorithm 
generalizes easily to sample from any distribution P{x) whose configuration space is of the 
form X"^. The algorithm can be described by the following pseudocode: 

1. Select Xq uniformly at random across A"". 

2. For i = 1 to t = T: 

3. Draw an integer i at random with uniform distribution across {1,2,..., n}. 



A, 



■rs 




^-l3{Es-Er) 



1 Q-l3{Es-Er 
Ps 



(416) 
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4. For each x e X, calculate 

P{X^ = = = ^ ''%rj'''7^^f \.y (417) 

5. Set a^t+i = a^t for all j i and = X*, where is drawn according to 

6. end 

7. Return the sequence Xf, t — 1,2, . . . ,T. 

It can be easily seen that the resulting Markov chain satisfies detailed balance and that in 
the case of binary alphabet (spin array) it implements the above expression of A^g- One can 
also easily generalize the Metropolis algorithm, in the same spirit, as e~^^^'~'^^^ is nothing 
but the ratio Ps/Pr- 
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